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Abstract 

Fractional kinetic equations of the diffusion, diffusion-advection, and Fokker-Planck type are presented 
as a useful approach for the description of transport dynamics in complex systems which are governed by 
anomalous diffusion and non-exponential relaxation patterns. These fractional equations are derived 
asymptotically from basic random walk models, and from a generalised master equation. Several physical 
consequences are discussed which are relevant to dynamical processes in complex systems. Methods of 
solution are introduced and for some special cases exact solutions are calculated. This report demonstrates 
that fractional equations have come of age as a complementary tool in the description of anomalous 
transport processes. © 2000 Elsevier Science B.V. All rights reserved. 

PACS: 05.40. - a; 05.40.Fb; 02.50.Ey 

Keywords: Anomalous diffusion; Fractional diffusion equation; Fractional Fokker-Planck equation; Anomalous relax- 
ation; Mittag-Leffler relaxation; Dynamics in complex systems 
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For the beginning is assuredly 

the end - since we know nothing, pure 

and simple, beyond 

our own complexities 

William Carlos Williams, Patterson 



1. Prologue: the scope, and why bother at all 

1.1. What can fractional equations do, what can they do better, and why should one care at all! 

Before we start off with the Introduction, we would like to address some points which we believe 
to strike many colleagues who are not familiar with the topic. 

The universality: The detailed structure of the propagator W(r, t), i.e., the probability density 
function (pdf) for the initial condition lim,^ + W(r,t) = 5(r), depends, in general, on the special 
shape of the underlying geometry. However, the interesting part of the propagator has the 
asymptotic behaviour log W(r, t) ~ — c£ u where £ = r/t a ' 2 > 1 which is expected to be universal. 
Here, w = 1/(1 — a/2) with the anomalous diffusion exponent a defined below. The fractional 
equations we consider in the following are universal in this respect as we do not consider any form 
of quenched disorder. Our results for anomalous diffusion are equivalent to findings from random 
walk models on an isotropic and homogeneous support. 

The non-universality: In contrast to Gaussian diffusion, fractional diffusion is non-universal in 
that it involves a parameter a which is the order of the fractional derivative. Obviously, nature 
often violates the Gaussian universality mirrored in experimental results which do not follow the 
Gaussian predictions. Fractional diffusion equations account for the typical "anomalous" features 
which are observed in many systems. 

The advantage to random walk models: Within the fractional approach it is possible to include 
external fields in a straightforward manner. Also the consideration of transport in the phase space 
spanned by both position and velocity coordinate is possible within the same approach. Moreover, 
the calculation of boundary value problems is analogous to the procedure for the corresponding 
standard equations. 

The comparison to other approaches: The fractional approach is in some sense equivalent 
to the generalised master equation approach. The advantage of the fractional model again lies in 
the straightforward way of including external force terms and of calculating boundary value 
problems. Conversely, generalised Langevin equations lead to a different description as they 
correspond to Fokker-Planck equations which are local in time and which contain time-dependent 
coefficients. In most cases of Brownian transport, the deterministic Fokker-Planck equation is 
employed for the description of stochastic dynamics in external fields. In analogy, we promote 
to use the fractional Fokker-Planck equation for situations where anomalous diffusion underlies 
the system. 

The mathematical advantage: A very convenient issue is that standard techniques for solving 
partial differential equations or for calculating related transport moments also apply to fractional 
equations which is demonstrated in the text. 
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The relation between the fractional solution and its Brownian counterpart There exists a trans- 
formation which maps the Brownian solution onto the corresponding fractional solution, an 
interesting relation which is useful for both analytic and numerical analysis. 

It is a simple approach: The appearance of fractional equations is very appealing due to their 
proximity to the analogous standard equations. It has been demonstrated recently that the 
fractional Fokker-Planck equation can be derived from a Langevin equation with Gaussian white 
noise for systems where trapping occurs. This offers some insight into the physical mechanisms 
leading to fractional kinetics. 

Fractional kinetic equations are not just another way of presenting old stories. We believe that 
they are a powerful framework which is of use for many systems. By relating our and others' work in 
that field and putting it in some more general context, the present report may be the basis for some 
active research on complex dynamics using a tool which is as old and new as classical calculus. 

1.2. What is the scope of this report! 

We do not present a report on anomalous transport theory: Anomalous diffusion is an involved field 
with intriguing subtleties. Accordingly, parameters and exponents can change in the course of time 
or when an external force is switched on or off, etc. In our considerations we do not touch on these 
issues but assume a given anomalous diffusion exponent. The present knowledge on anomalous 
diffusion is compiled in a number of review articles quoted in the list of references. 

The focus is on subdijfusion: The major body of the text is concerned with transport processes 
which are, in the force-free limit, slower than Brownian diffusion. In a lose sense we also call such 
processes subdiffusive or dispersive which take place in an external force field and whose force-free 
limit corresponds to subdiffusion. Occasionally, such processes are called fractional. 

Some care be taken with Levy flights: Levy flights do not possess a finite mean squared displacement. 
Their physical significance therefore has been questioned as particles with a finite mass should not 
execute long jumps instantaneously. We do consider them to some extent in the following as there are 
special cases whose description in terms of Levy flights corresponds to physical principles. 

Levy walks: The proper way of considering systems which feature Levy type jump length 
distributions is to introduce a finite speed for the test particle, a model referred to as Levy walk. 
A first step towards a fractional dynamics formulation has been presented recently. 

Since Newton integer order calculus has been used in physical modelling. Newton's rival Leibniz' 
prophesy already claimed in 1695: "Thus it follows that d 1/2 x will be equal to x^/dx : x, an apparent 
paradox, from which one day useful consequences will be drawn" - we believe that the following 
account represents an affirmative answer to this first statement towards a fractional calculus. 



2. Introduction 

2.1. Anomalous dynamics in complex systems 

Complex systems and the investigation of their structural and dynamical properties have 
established on the physics agenda. These "structures with variations" [1] are characterised through 
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(i) a large diversity of elementary units, (ii) strong interactions between the units, or (iii) a non- 
predictable or anomalous evolution in the course of time [2]. Complex systems and their study 
play a dominant role in exact and life sciences, embracing a richness of systems such as glasses, 
liquid crystals, polymers, proteins, biopolymers, organisms or even ecosystems. In general, the 
temporal evolution of, and within, such systems deviates from the corresponding standard laws 
[3,4]. With the development of higher experimental resolutions, or the combination of different 
experimental techniques, these deviations have become more prominent, as the accessible larger 
data windows are more conclusive. 

Thus, relaxation in complex systems deviates from the classical exponential Debye pattern [5-10] 

&(t) = O exp( - t/x) , (1) 

and can often be described in terms of a Kohlrausch-Williams-Watts stretched exponential law 
[11-15] 

&(t) = O exp( - (t/xf) (2) 
for < a < 1, or by an asymptotic power law [9,13,14,16-19] 

0(t) = <2> (1 + f/T) - " (3) 

with n > 0. It is also possible to observe transitions from the stretched exponential pattern (2) to the 
power-law behaviour (3) [20-22]. 

Similarly, diffusion processes in various complex systems usually no longer follow Gaussian 
statistics, and thus Fick's second law fails to describe the related transport behaviour. Especially, 
one observes deviations from the linear time dependence of the mean squared displacement 

<x 2 (t))~K lt (4) 

which is characteristic of Brownian motion, and as such a direct consequence of the central limit 
theorem and the Markovian nature of the underlying stochastic process [23-39]. Instead, anomal- 
ous diffusion is found in a wide diversity of systems, its hallmark being the non-linear growth of the 
mean squared displacement in the course of time. In this report, we concentrate on the power-law 
pattern 

<x 2 (0> ~ K a f (5) 

which is ubiquitous to a diverse number of systems [14,38-55]. There exists a variety of other 
patterns such as a logarithmic time dependence which we do not touch upon here. The anomalous 
diffusion behaviour manifest in Eq. (5) is intimately connected with the breakdown of the central 
limit theorem, caused by either broad distributions or long-range correlations. Instead, anomalous 
diffusion rests on the validity of the Levy-Gnedenko generalised central limit theorem for 
such situations where not all moments of the underlying elementary transport events exist 
[38,42,47,56-58]. 1 Thus, broad spatial jump or waiting time distributions lead to non-Gaussian 



1 Actually, the notion of diverging moments goes back to the formulation of the "St. Petersburg paradox" by Nicolaas 
Bernoulli and its analysis by Daniel Bernoulli in 1724 [38,49]. 
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Fig. 1. Different domains of anomalous diffusion, defined through the mean squared displacement, Eq. (5), parametrised 
by the anomalous diffusion exponent a: (a) subdiffusion for < oc < 1, (b) superdiffusion for a > 1. On the threshold 
between sub- and superdiffusion is the normal Brownian diffusion located. Another special case is ballistic motion 
(a = 2). 



propagators and a possibly non-Markovian time evolution of the system, the latter being a typical 
manifestation of non-local temporal phenomena encountered in a broad range of systems 
[1,3,4,59-63]. Note that the generalised diffusion coefficient K a in Eq. (5) has the dimension 
= cm 2 s~*. According to the value of the anomalous diffusion exponent a, defined in Eq. (5), 
one usually distinguishes several domains of anomalous transport, as sketched in Fig. 1. In what 
follows, the main emphasis will be laid on the description of subdiffusive phenomena, which 
correspond to < a < 1. 

2.2. Historical remarks 

The stochastic formulation of transport phenomena in terms of a random walk process, 2 as well 
as the description through the deterministic diffusion equation 3 are the two fundamental concepts 
in the theory of both normal and anomalous diffusion. Indeed, the history of this dual description 
basing on erratic motion and on a differential equation for the probability density function is quite 
interesting and much worth a short digression. 

Thus, small flickering of coal dust particles on the surface of alcohol was observed by the Dutch 
physician Jan Ingenhousz 4 as early as in 1785. In 1827, the Scottish botanist Robert Brown [67] 
observed similar irregular movement of small pollen grain under a microscope. 5 At about the same 
time, in 1822, Joseph Fourier came up with the heat conduction equation, on the basis of which 



2 The random walk concept was formally introduced by Karl Pearson in 1905 [65]. 

3 We use the word deterministic in order to distinguish the partial differential equation for the probability density 
function W(x, t) from a stochastic differential equation like the Langevin equation [37,64], compare, e.g., Risken [36]. 

4 Later the personal practitioner of Austrian empress Maria Theresia, he also discovered the respiration process of 
plants. 

5 He also recognised the significance of the cell nucleus. 
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Fig. 2. Recorded random walk trajectories by Jean Baptiste Perrin [72], Left part: three designs obtained by tracing 
a small grain of putty {mastic, used for varnish) at intervals of 30 s. One of the patterns contains 50 single points. Right 
part: the starting point of each motion event is shifted to the origin. The figure illustrates the pdf of the travelled distance 
r to be in the interval (r, r + dr), according to (2ti£ 2 )~ 1 exp( — r 2 /[2£ 2 ])2nr dr, in two dimensions, with the length variance 
£ 2 . These figures constitute part of the measurement of Perrin, Dabrowski and Chaudesaigues leading to the determina- 
tion of the Avogadro number. The result given by Perrin is 70.5 x 10 22 . The remarkable ceuvre of Perrin discusses all 
possibilities of obtaining the Avogadro number known at that time. Concerning the trajectories displayed in the left part 
of this figure, Perrin makes an interesting statement: "Si, en effet, on faisait des pointes de seconde en seconde, chacun de 
ces segments rectilignes se trouverait remplace par un contour polygonal de 30 cotes relativement aussi complique que le 
dessin ici reproduit, et ainsi de suite". [If, veritably, one took the position from second to second, each of these rectilinear 
segments would be replaced by a polygonal contour of 30 edges, each itself being as complicated as the reproduced 
design, and so forth.] This already anticipates Levy's cognisance of the self-similar nature, see footnote 9, as well as of the 
non-differentiability recognised by N. Wiener. 



A. Fick set up the diffusion equation in 1855 [68]. 6 Subsequently, the detailed experiments by Gouy 
proved the kinetic theory explanation given by C. Weiner in 1863. After attempts of finding 
a stochastic footing like the collision model by von Nageli and John William Strutt, Lord 
Rayleigh's results, it was Albert Einstein who, in 1905, unified the two approaches in his treatises on 
the Brownian motion, a name coined by Einstein although he reportedly did not have access to 
Brown's original work. Note that a similar description of diffusion was presented by the French 



6 In the historical context, note that the theory of the continuum formulation of fluid dynamics had already been fully 
developed at that time. Thus, some of its milestones date back to the 18th and first half of the 19th century, such as 
Bernoulli's equation (1738); Euler's equation (1755); Navier's (1827) use as a phenomenological model and Stokes' (1845) 
derivation of the Navier-Stokes equation. Maxwells' dynamical theory of gases dates back to 1867 and Boltzmann's 
transport equation was published in 1872 for the description of collision processes. The latter is the footing for the 
atomistic random walk approach to Brownian motion. 
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mathematician Louis Bachelier in his 1900 thesis [69], in terms of stock values instead of physical 
quantities [38,70]. An important application of Einstein's results was the independent measure- 
ment of the Avogadro number 7 by Jean Baptiste Perrin, A. Westgren and Eugen Kappler 
[26,28,71-74], to a rather high accuracy. Some results of Perrin are shown in Fig. 2 and they are 
part of the work which won him the Nobel Prize in 1926. The random walk which can be 
experimentally observed, represents therefore a link between the microscopic dynamics of small 
atoms bombarding a larger particle in suspension, and macroscopic observables like the diffusion 
coefficient, or the Avogadro number. In Fig. 3, we reproduced data obtained by Kappler with his 
high-accuracy set-up using an optical detection method (a detailed explanation is given in the 
figure caption). Einstein's ideas also set the scene for Langevin's treatment [37,64] of Brownian 
motion with the assumption of an external erratic force, and the Fokker-Planck [78,79], 
Smoluchowski 8 and Klein-Kramers [81,82] theories which culminated in the treatises of Ornstein 
and Uhlenbeck, Chandrasekhar and others, and later in the works of Elliott Montroll, and 
collaborators [24,83-86]. 

The mathematical treatment of Brownian motion is mainly due to Norbert Wiener who proved 
that the trajectory of a Brownian particle is (almost) everywhere continuous but nowhere differenti- 
able [87]. This observation is related to the self-affine nature of the diffusion process whose 
resulting spatial trajectory is self-similar 9 [27,38,88-93]. Further important mathematical contri- 
butions attribute to J. L. Doob, Mark Kac, W. Feller, and others. 

2.3. Anomalous diffusion: experiments and models 

Anomalous diffusion has been known since Richardson's treatise on turbulent diffusion 
in 1926 [94]. Within transport theory it has been studied since the late 1960s. In particular, 
its theoretical investigation was instigated by Scher and Montroll in their description of dispersive 
transport in amorphous semiconductors, a system where the traditional methods proved to 
fail. The predictions of their continuous time random walk approach were very distinct from 
its Brownian counterpart and were shown to provide explanations for a variety of physical 
quantities and phenomena in numerous experimental realisations [95-99]. Important contribu- 
tions are also due to Weiss [39] and Shlesinger [100]. Besides the random walk description, 
generalisations of the diffusion equation were developed which account for the anomalous 
transport statistics. 

Today, the list of systems displaying anomalous dynamical behaviour is quite extensive 
[14,40-55,101,102] and hosts, among others, the following systems in the subdiffusive regime: 



7 Historically, often referred to as the Loschmidt number [28]. 

8 The monovariate Fokker-Planck equation discussed here is often referred to as Smoluchowski equation [80], 
9 "Le processus stochastique, que nous appellerons mouvement brownien lineaire, est une schematisation, qui 

represente bien les proprietes du mouvement brownien reel observables a une echelle assez petite, mais non infiniment 
petite, et qui suppose que les memes proprietes existent a n'importe aquelle echelle". (Paul Levy [27]) [The stochastic 
process which we will call linear Brownian motion is a schematic representation which describes well the properties of 
real Brownian motion, observable on a small enough, but not infinitely small scale, and which assumes that the same 
properties exist on whatever scale.] 
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charge carrier transport in amorphous semiconductors [95-99,101,103,104], nuclear magnetic 
resonance (NMR) diffusometry in percolative [105,106], and porous systems [107,108], Rouse or 
reptation dynamics in polymeric systems [109-115], transport on fractal geometries [116,117], the 
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Fig. 3. Erroneous behaviour of Brownian motion made visible in a high-precision measurement. Data from an 
Edelmann recorder, obtained by Kappler in 1931 [74]. Following calculations of Smoluchowski [75] and first 
measurements by Gerlach and Lehrer [76,77], Kappler monitored the Brownian motion of a small mirror (surface 
approx. 1mm 2 ), suspended from a fine quartz thread (several cm long and some um thick). The mean square of 
the torsional displacement, cp 2 , follows the relation Dcp 2 = k B T, where D is the directional force of the suspension [74]. 
The facsimiles show four different realisations. From his data, Kappler obtained the Avogadro-Loschmidt number 
N L = 60.59 x 10 22 + 1%, to a remarkable accuracy. 



R. Metzler, J. Klafter / Physics Reports 339 (2000) 1-77 



11 



Registrieraufiialime der Brownachea Bewegang (natiirltclie Grofie). 
Direktionskraft 9,428 • 10 — 6 abs. Einh, Tragheitamoment: 1-10~* aba. Einh. Abatand Spiegel-Kamera : 72,1 cm. 
Zeitmarke: 30 sec dx = 1 mm. a) Atmoaphiirendruck. Temperatur 13" C 

Fig. 5 a 



r A. .a ' • 



liegistrieraufnahme der Brownselien Bewc^ung (naturliclic GrSBe). 
Dircktioiiakraft 9,428- 10 -3 aba. Einh. Triigheitsmoment 1 • 10^ 7 aba. Einh. Abstand Hptegel-Kamcra: 72,1cm. 
Zeitmarke: 30 aec dx = 1 mm. b) 1-10 -3 mm Hg. Temperatur 13° C 
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Fig. 3 (continued). 



diffusion of a scalar tracer in an array of convection rolls [118,119], or the dynamics of a bead in 
a polymeric network [120,121]. Superdiffusion or Levy statistics are observed in special domains of 
rotating flows [122], in collective slip diffusion on solid surfaces [123], in layered velocity fields 
[124,125], in Richardson turbulent diffusion [94,126-129], in bulk-surface exchange controlled 
dynamics in porous glasses [130-132], in the transport in micelle systems and in heterogeneous 
rocks [133-135], in quantum optics [136,137], single molecule spectroscopy [138,139], in the 
transport in turbulent plasma [140], bacterial motion [141-145] and even for the flight of an 
albatross [146]. 

Anomalous diffusion in the presence or absence of an external velocity or force field has been 
modelled in numerous ways, including (i) fractional Brownian motion dating back to Benoit 
Mandelbrot [89-93,147], (ii) generalised diffusion equations [148], (iii) continuous time random 
walk models [52,95-99,101,102,149-155], (iv) Langevin equations [156-160], (v) generalised Lan- 
gevin equations [62,161-163], (vi) generalised master equations [164-167], or (vii) generalised 
thermostatistics [168-172]. For anomalous diffusion, only the approaches (iii) and (v) incorporate 
the system's memory and the special form which is to be expected for the pdf, in a consistent way. 
The disadvantage in the continuous time random walk and the generalised master equation 
approaches is that there is no straightforward way to incorporate force fields, boundary value 
problems, or to consider the dynamics in phase space. 

The alternative approach to anomalous kinetics which we are going to present is given in 
terms of fractional equations which appear to be tailored for such kind of problems like the 
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consideration of external fields and boundary value problems. In the original works [173-175] 10 it 
was realised that the replacement of the local time derivative in the diffusion equation by 
a fractional operator accounts for the memory effects which are connected with many complex 
systems. 

Recently, a decade after their introduction, such fractional kinetic equations have attracted much 
interest. They are presently being extensively studied and recognised as important tools in the 
description of anomalous transport processes in both absence and presence of external velocity or 
force fields. Especially in the latter case, their mathematical structure allows for the application of 
known methods of solution. In the course of this development, a number of works has been 
published dealing with fractional relaxation equations and fractional rheological models 
[20,21,176-182], fractional diffusion equations (FDEs) [173-175,183-201], fractional diffusion- 
advection equations (FDAEs) [202-209], and fractional Fokker-Planck equations (FFPEs) 
[158-160,193,202,203,210-221]. Thereby, various generalisations to fractional order have been 
employed, i.e. different fractional operators have been introduced to replace either the time 
derivative or the occurring spatial derivatives, or both. 

In first attempts of generalising the standard diffusion equation for the description of diffusion 
processes on fractal geometries of dimension d f , radius-dependent diffusion coefficients were 
assumed. O'Shaugnessy and Procaccia studied the generalised diffusion equation [148] 11 

d ^ = R 1 - d <^-R i ' - 'K(R)^-W(R, t) (6) 
dt dR 'dR ' ' 

with the radius-dependent diffusion coefficient K(R) = KR , and derived the corresponding 
propagator 

W(R,t) = A(0,d f )(K\_2 + 0] 2 O-^ /(2 + 0) exp( - K ^° Q)2 ^ • (?) 

Here, = d f + a — 2 is connected to the power-law index a of the radius-dependent, integrated 
electrical resistance, {%(R) ~ R a of the underlying fractal structure the mass density of which scales 
as oc R df ~ 3 in a 3-dimensional embedding. The mean squared displacement for this process, given 
by Eq. (5) with a = 2/(2 + 0), can readily be inferred. As is positive, this result implies 
subdiffusion [148]. Further investigations indicated that the asymptotic form of the propagator on 
fractals such as the Sierpihski gasket, is given by the scaling form [222,223] 

W(R,t)^A(aJ)^cxp(-c^) (8) 

with £ = R/t* 12 , u = 1/(1 — a/2), and /? is a system-dependent quantity [224]. Eq. (8) is at variance 
with the above result, Eq. (7). For subdiffusion, < a < 1, which prevails on fractal structures, 
we (0,2) so that expression (8) is often referred to as a stretched Gaussian. Conversely, fractional 



10 In Ref. [173], note that the author introduces the defining expression for a fractal operator, but does not make 
explicit use of fractional calculus, neither of the term fractional as such. 

1 1 The capital R refers to the radius averaging over the fractal support which is necessary in order to obtain the smooth 
radius dependence. 
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diffusion equations for transport on fractal structures were shown to comply with the basic 
properties, like the propagator (8), the returning probability to the origin, the mean squared 
displacement, and the non-Markovian nature [183,184,225]. 

For homogeneous, isotropic systems which we are interested in, one knows from random walk 
models that the propagator behaves asymptotically like 

W(r, t) ~ b r ydl2 ^ exp( - b t <T) (9) 

with £ = d ll2 r/t a/2 , v = 1/(1 — a/2), /? = — d(l — a)/(2 — a), and b and b 1 are constants depending 
on a and the dimension d [238,239]. Such behaviour is reproduced by the fractional diffusion 
equations that were anticipated by Balakrishnan [173] and first formulated by Wyss [174], and 
Schneider and Wyss [175]. In Section 3 we show that the asymptotic behaviour (9) is consistent 
with the asymptotic expansion of the exact solution of the fractional diffusion equation. Fractional 
kinetic equations, their foundation and solution form the centre-piece of this report. 

In what follows, we present the basic principles and physical properties connected with fractional 
kinetic equations. We show that fractional equations arise naturally in the diffusion limit of certain 
random walk schemes. By discussing methods of solution and deriving explicit solutions, we 
demonstrate the usefulness of the fractional approach. At first, we introduce FDEs basing on the 
continuous time random walk model where the transport events are subject to broad statistics. 
Subsequently, FDAEs and FFPEs are presented and they describe the transport in an external 
velocity or force field. In the final section, a physical derivation on the basis of a Langevin equation 
with white Gaussian noise is discussed which leads to a fractional Klein-Kramers equation. From 
the latter, the FFPE is consistently recovered. In the Appendices we have compiled some basic 
definitions and useful relations which are of relevance in the main text. Thus, we give a primer on 
fractional calculus and the special functions which emerge when dealing with fractional differential 
equations, as well as a short introduction to Levy stable laws. Finally, we list the abbreviations 
and the notation used. Note that throughout the text we denote the Laplace and Fourier 
transforms of a function by stating the explicit dependence on the associated variable, e.g. 
W(k,u) = 3?{£e{W{x,t);t^u};x^>k}. 

The numerous illustrations spread throughout the text are to visualise the often striking 
differences in functional behaviour of the normal and anomalous cases, especially the perseverance 
of the initial condition in the subdiffusive domain. 

The discussion in the remainder of this report is restricted to the one-dimensional case, with special 
emphasis on subdiffusion phenomena. The equations which describe subdiffusion presented in the 
text can be extended to higher dimensions through a replacement of the derivatives in respect of the 
position coordinate by corresponding orders of the V operator. Some remarks on higher-dimensional 
systems are contained in Ref. [215]. For those equations which describe situations with Levy 
distributed jump lengths and which therefore contain a generalised Laplacian, we refer to the 
definitions in Ref. [226] and the discussions in Refs. [156-158] for the multi-dimensional case. 



3. From continuous time random walk to fractional diffusion equations 

In our quest of establishing the fractional diffusion equation (FDE), we borrow from the ideas of 
connecting the random walk approach with the continuum description through the diffusion 
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Fig. 4. Schematic representation ol a Brownian random walk. The walker jumps at each time step f = 0, 
At, 2At, . . . , n At, . . . to a randomly selected direction, thereby covering the distance Ax, the lattice constant. 



equation, and we start off with the continuous time random walk scheme. The latter is flexible 
enough to account for the rich panel of transport regimes encountered in complex systems. After 
establishing the fundamental framework of random walks and recovering the standard diffusion 
equation, we move on to the discussion of the continuous time random walk framework, and the 
derivation of the FDE. This equation will be shown to enable an investigation of subdiffusive 
phenomena, and Levy flights with the tools well-known from dealing with the standard diffusion 
equation. 

3.1. Revisiting the realm of Brownian motion 

A typical Brownian walk like Perrin's original data, is schematically displayed on a two- 
dimensional lattice in Fig. 4. In discrete time steps of span At, the test particle is assumed to jump to 
one of its nearest neighbour sites, here displayed on a square lattice with lattice constant Ax, the 
direction being random. Such a process can be modelled by the master equation 

Wj(t + At) = Wj-At) + W J+ i(t) (10) 

in the one-dimensional analogue, as the process is local in both space and time. In Eq. (10), the 
index denotes the position on the underlying one-dimensional lattice. Eq. (10) defines the pdf to be 
at position j at time t + At in dependence of the population of the two adjacent sites j ± 1 at time t. 
The prefactor 1/2 expresses the direction isotropy of the jumps. In the continuum limit At -> and 
Ax -> 0, Taylor expansions in At and Ax, 



Wj(t + At) = Wj(t) + At— + + 0([At] 2 ) 



(11) 
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and 



c)W (Ax) 2 c) z W 

Wj±1 (t) = W(x,t) ± Ax— + L_L _ + Q([Ax] 3 ) , (12) 



lead to the diffusion equation 

QW _ Q 2 
~8T ~~ J 8x 



Ki— 2 W(x,t), (13) 



on taking along the lowest orders in At and Ax. The continuum limit thereby has to be drawn such 
that the quotient 

K, = lim ^ (14) 

is finite. K 1 is called the diffusion constant and is of dimension [Ki] = cm 2 s _1 . 

The diffusion equation (13) is one of the most fundamental equations in physics, being a direct 
consequence of the central limit theorem [27,42,58]. Under the condition that the first two 
moments of the pdf, describing the appropriately normalised distance covered in a jump event and 
the variance, X = and X 2 , as well as the mean time span At between any two individual jump 
events, exist, the central limit theorem assures that the random walk process is characterised by 
a mean velocity V = X/At and a diffusion coefficient K = (2At) _1 [X 2 - X 2 ] [27,42]. Further- 
more, for long times, i.e., a large enough number of steps, the pdf of being at a certain position x at 
time t, is governed by the diffusion equation (13), and it is given by the Gaussian shape 

WM= ^hrA-^)- ,15) 

W(x, t) from Eq. (15) is called the propagator, i.e., the solution of the diffusion equation (13) for the 
sharp initial condition W (x) = lim ( ^ 0+ W(x, t) = S(x). Individual modes of Eq. (13) decay expo- 
nentially in time, 

W{k,t) = Qxv{-K Y k 2 t) , (16) 
with the Fourier transformed diffusion equation, 
dW 

— = -K t k 2 W(k,t), (17) 



being a relaxation equation, for a fixed wavenumber k. 



3.2. The continuous time random walk model 



For the generalisations to anomalous transport, we choose the continuous time random walk 
(CTRW) scheme as the starting point. In parallel to the complementary, dual approach in the 
standard diffusion problem, we will then develop a generalised diffusion equation of fractional 
order on the basis of the CTRW. 
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Fig. 5. Continuous time random walk (CTRW) model. Left: CTRW process on a two-dimensional lattice, generalising 
the Brownian situation from Fig. 4. The waiting times are symbolised by the waiting circles the diameter ol each is 
proportional to the waiting time which is to be spent on a given site before the next jump event occurs. The jump lengths 
are still equidistant. Right: (x, t) diagram of a one-dimensional CTRW process where both jump lengths and waiting times 
are drawn from pdfs which allow for a broad variation of the corresponding random variables. 



The CTRW model is based on the idea that the length of a given jump, as well as the waiting time 
elapsing between two successive jumps are drawn from a pdf \j/(x, t) which will be referred to as the 
jump pdf. From \j/(x, t), the jump length pdf 



;.(x) = J dt\Kx,t) (18) 
and the waiting time pdf 



w(t) = 



dx ij/(x, t) (19) 



can be deduced. Thus, X{x)dx produces the probability for a jump length in the interval (x, x + dx) 
and w(t)dt the probability for a waiting time in the interval (t, t + dt). If the jump length and 
waiting time are independent random variables, one finds the decoupled form ij/(x, t) = w(t)X(x) for 
the jump pdf \j/(x, t). If both are coupled, i.e., \j/(x, t) = p(x\t)w(t) or \f/(x, t) = p(t\x)A(x), a jump of 
a certain length involves a time cost or, vice versa; i.e., in a given time span the walker can only 
travel a maximum distance. In what follows, we employ the decoupled version. A schematic 
cartoon of the CTRW model is drawn in Fig. 5. 

Different types of CTRW processes can be categorised by the characteristic waiting time 



T = J dt w(t)t (20) 
and the jump length variance 

'OO 

I 2 = dxX{x)x z (21) 
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being finite or diverging, respectively. With these definitions, a CTRW process can be described 
through an appropriate generalised master equation [150,151,164,166,167], via a set of Langevin 
equations [159,160,216,217], or by the equation [151] 



n(x, t) = dx' dt' tj(x', t')\j/(x -x',t - t') + 8(x)8(t) (22) 

J - qo Jo 

which relates the pdf rj(x, t) of just having arrived at position x at time t, with the event of having 
just arrived at x' at time t', t](x', t'). The second summand in Eq. (22) denotes the initial condition of 
the random walk, here chosen to be 8(x). Consequently, the pdf W(x, t) of being in x at time t is 
given by 



W(x,t) 



dtr]{x, t')W{t - t') , (23) 



i.e., of arrival on that site at time f, and not having moved since. The latter is being defined by the 
cumulative probability 

«P(t) = 1 - | dt' w(t') (24) 

assigned to the probability of no jump event during the time interval (0, t). In Fourier-Laplace 
space, the pdf W(x, t) obeys the algebraic relation [151] 

mu) = '-("> ■"•<*> , (25) 

u 1 — yj(k, u) 

where W (k) denotes the Fourier transform of the initial condition W (x). 



3.3. Back to Brownian motion 

Consider now different cases of the CTRW model defined through the decoupled jump pdf 
\j/{x, t) — w(t)X(x). If both characteristic waiting time and jump length variance, T and I 2 , are finite, 
the long-time limit corresponds to Brownian motion. Let us consider, for instance, a Poissonian 
waiting time pdf w(t) = T~ 1 exp( — t/x) with T = t, together with a Gaussian jump length pdf 
X(x) = (47t<T 2 )~ 1/2 exp( — x 2 /(4a 2 )) leading to I 2 = 2a 2 . Then, the corresponding Laplace and 
Fourier transforms are of the forms 

w(u) ~ 1 - m + 0(t 2 ) , (26) 

X(k) ~ 1 - a 2 k 2 + 0(ic 4 ) . (27) 

In fact, any pair of pdfs leading to finite T and I 2 leads to the same result, to lowest orders, and 
thus in the long-time limit [151]. Installing Eqs. (26) and (27) into Eq. (25), one readily recovers, for 



18 



R. Metzler, J. Klafter / Physics Reports 339 (2000) 1-77 



the initial condition W (x) = 8(x), the Fourier- Laplace space transform of the propagator. 
1 

u + K±k 



W(k,u) = - , - ,, 2 , (28) 



with K 1 = (t 2 /t. Back-transformed to (x, ^-coordinates, this is but the well-known Gaussian 
propagator, Eq. (15). After multiplication with the denominator in Eq. (28), and making use of 
the differentiation theorems of Fourier (i.e., ^{d z W(x, t)/dx 2 } = — k 2 W(k,t)) and Laplace 
(i.e., J?{dW(x, t)/dt} = uW(x,u) — W (x)) transformations, Fick's second law (13) is immediately 
obtained. Note that the notion long-time, equivalent to the diffusion limit, is only relative in respect 
to the time scale t. In Fourier-Laplace space, the diffusion limit is given through the assumption of 
(0,0) [150,223,227,228]. 

3. 4. Long rests: a fractional diffusion equation describing subdiffusion 

Consider the following situation, sometimes referred to as fractal time random walk [101], where 
the characteristic waiting time T diverges, but the jump length variance I 2 is still kept finite. To 
this end, a long-tailed waiting time pdf with the asymptotic behaviour [83-85] 

w(t) ~ A x (T/t) 1+ « , (29) 

for < a < 1 is introduced, which has the corresponding Laplace space asymptotics 
[83-85,151,229-230] 12 

w(u) ~ 1 - (mf . (30) 

Again, the specific form of w(t) is of minor importance. Consequently, together with the Gaussian 
jump length pdf characterised through Eq. (27), the pdf in Fourier Laplace space becomes 



\W (k)lu\ 
1 + K x u-*k 2 



W(k,u) = ,\ (31) 



in the (k, u) -> (0, 0) diffusion limit. Employing the integration rule for fractional integrals [226,232], 
^{ D t -"W(x,t)} = u- p W(x,u), p > , (32) 
one infers the fractional integral equation 

W(x, t) - W (x) = oDrK x ^ W(x, i) (33) 

from relation (31). By application of the differential operator d/dt, one finally arrives at the FDE 



QW 



^2 



8, oA'-K^WM. (34) 



12 Note that \jj(u = 0), i.e. lim„^ jo dte u '\j/(t) is but the normalisation of the waiting time pdf, i.e. \j/{u = 0) = 1. 
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The Riemann-Liouville operator D} a = (d/dt) D t a , for < a < 1, is defined through the 
relation [226,232-235] 



1 8 

7 in.v./i 



V, W(x, f) 
o (t-0 



* 7^9£h • (35) 



r(a)dt 

Its fundamental property is the fractional order differentiation of a power, 

o^-t^S^t'*- 1 . (36) 
r(p + a) 

In fact, it can be shown that the more general relation 

holds, for any real q. Especially, the Riemann-Liouville fractional differintegration of a constant 
becomes 

The special cases of integer order differentiation of a constant, d"l/df" = 0, are included through 
the poles of the Gamma function for q = 1,2, ... .A more detailed introduction to fractional 
differintegration is given in Appendix A. 

Thus, the integrodifferential nature of the Riemann-Liouville fractional operator D} ~* accord- 
ing to Eq. (35), with the integral kernel M(t) oct 1-1 , ensures the non-Markovian nature of the 
subdiffusive process defined by the FDE (34). Indeed, calculating the mean squared displacement 
from relation (31) via the relation <x 2 > = lim k ^ { — (d 2 /dk 2 )W(k,u)} and subsequent Laplace 
inversion, the result 

2K 
T(l + a) 

is obtained. Alternatively, it can be inferred from the FDE (34) through integration over J^^dx x 2 , 
leading to (d/dt)<x 2 (t)> = £> t 1_a 2K a = 2X a f" _1 /r(a). 
Rewriting the FDE (34) in the equivalent form 

D?W - jJ^Woix) = K a ^W(x,t) , (40) 

the initial value W (x) is seen to decay with the inverse power-law form (t -a /-T(l — ct))W (x), and 
not exponentially fast as for standard diffusion [215]. Note that in the limit a -> 1, the FDE (34) 
reduces to Fick's second law, as it should. The generalised diffusion constant K a which appears in 
the FDE (34), is defined by 



K x = a 2 lx 



(41) 
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in terms of the scales <x and t, leading to the dimension [K a ] = cm s~*. The FDE (34) was first 
considered in the integral form (33) by Schneider and Wyss [175]. An equivalent form was 
considered by Balakrishnan [173], and a differential form by Wyss [174]. 

A closed-form solution for the FDE (34) can be found in terms of Fox functions, the result being 



W(x, t) = 



1 



J2l,2 



" x 2 


(1 — a/2, a) 


|_4K a f" 


(o,i), (U) _ 



(42) 



introducing the Fox function H\ ;°. The Fox function is defined in Appendix B. Note that the result 
(42) can be rewritten in the alternative form 



W(x, t) = 



1 



(1 - a/2, a/2)" 
(0,1) 



(43) 



employing the definition of the Fox function and the duplication rule of the Gamma function 
[236]. 

Due to the occurrence of non-integer powers of the Laplace variable u in the expression for 
W(k, u), Eq. (31), a direct Laplace inversion is not tabled. There are three basic methods to compute 
the inversion: (i) First applied by Wyss [174], and Schneider and Wyss [175], the Mellin technique 
can overcome this problem by the roundabout way through Mellin space. Thereby, the path 
integral defining the Mellin inversion has a similar structure as the definition of the Fox functions, 
Eq. (B.8), so that the result can be directly inferred from its Mellin transform, (ii) One can identify 
the expression for W(k,u), Eq. (31), with its corresponding Fox function, and then use the existing 
rules for the Fox functions to calculate the Laplace and Fourier inversions, see Refs. [180,237]. 
The result is again a Fox function, which can be simplified by standard rules, to obtain the above 
results, (iii) One can first Fourier invert W(k, u), to obtain 



W(x, u) = \u al2 -^ exp( - \x\u al2 ) 



(44) 



expand the exponential function in its Taylor series, and invert term-by-term, using the rule (37). 
The final result is a power series, which can be shown to be identical with expression (43) [215]. 
Without the identification as a Fox function, the obtained series does not render any straightfor- 
ward information on the stretched exponential asymptotics (45) derived below from standard 
properties of the Fox function. 

Employing some standard theorems of the Fox function, one can derive the asymptotic stretched 
Gaussian behaviour 



W(x,t) 



1 



1 / 2 \a-«)/(2-«y \ x \ 



4nKj*\l 2 - a\a 

2 - a/a^ 2 "** 



■(l-«)/(2-«) 



x exp 



2 \2 



K a t« 

1/(1 -a/2) 



'K„f 



(45) 



valid for |x| > ^/K x t a . The functional form of the result (45) is equivalent to the CTRW findings 
reported by Zumofen and Klafter [238,239], see Eq. (9). Furthermore, W(x, t) can be represented 
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Fig. 6. Propagator W(x,t) for subdiffusion with anomalous diffusion exponent a = 1/2, drawn for the consecutive times 
t = 0.1, 1, 10. The cusp shape of the pdf is distinct. 



via the series expansion 
W(x, t) -- 1 



I 



1)" 



4K a fn%n\r(l-<x[n + l]/2)\K a f 



nil 



(46) 



in computable form. If a is a rational rather than a real number, the Fox function in Eqs. (42) and 
(43) can be simplified. Thus, for a = 1/2, it can be rewritten in terms of the Meijer G-function 



W(x, t) 



1 



2n 2 K 1/2 t 
1 



1/2 



1/2' 



(-1/2 



1/2' 



.1/2 



ej2,0 
-"0,2 



t/3,0 
^0,3 



/"3,0 



8K 1/2 ^ 2 

v 2 



(0, i),(ii) 



(0,1) &!),&!) 



16K 



1/2' 



a/2 



n 1 1 

u ? 4? 2 



(47) 



by twice using the duplication formula of the Gamma function [236] in the Mellin-Barnes type 
integral, Eq. (B.8), defining the Fox function [237,240-243]. This representation is useful, as 
the Meijer G-function belongs to the implemented special functions of Mathematica [244], the 
notation being 



1 



7^0,3 



x 2 \ 2 

16^] 



nii 

U> 4> 2 



= l/(8Pi~3tf(l/2)> 

MeijerG[{{}, {}}, {{0,1/4,1/2}, {}},x~4/(16~2 t)] 



(48) 
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Fig. 7. Propagator W(x, t) for Brownian diffusion (a = 1) for the times t = 0.05, 0.2 and 1. The much smoother shape 
close to the origin sets this solution apart from the subdiffusive counterpart drawn in Fig. 6. 



That way, the graphical representation of W(x, t) for the subdiffusive case a = 1/2 is obtained, 
which is displayed in Fig. 6. In comparison to the standard Gaussian result, shown in Fig. 7, the 
pronounced cusps of the subdiffusive propagator are distinct. Note that single modes of the FDE 
(34) decay in accordance to the Mittag-Leffler pattern 



with the asymptotic power-law behaviour W(k, t) ~ (K x k 2 t a r(l — a)) -1 . This typical Mittag- 
Leffler behaviour of the mode relaxation replaces the exponential mode relaxation (16) occurring 
for normal diffusion, and it is discussed in more detail in Section 5. The Mittag-Leffler function is 
introduced in Appendix B. 

Recently, it has been shown how boundary value problems for the FDE (34) can be solved [200]. 
As the jump length pdf in the subdiffusive case < a < 1 is narrow, i.e., it possesses a finite variance 
I 2 , one can apply the method of images due to Lord Kelvin which is summarised in the book of 
Feller [88]. Consider, for instance, the half space problem with a reflecting boundary at the origin. 
This situation is defined through the von Neumann condition, (dQ(x, t)/dx)\ x = = 0, where Q, 
specified below, denotes the image solution of this boundary value problem. Suppose the initial 
condition to be a sharp distribution at x , W (x) = 8(x — x ). Then the free solution can be 'folded' 
along a line through the origin, perpendicular to the x axis, i.e., the unrestricted solution is taken, 
and the portion which spreads to the space region opposite to x , in respect to the origin, is 
reflected at this line; the final result fulfils the von Neumann condition. The solution is thus given 
by the function [88] 



W(k,t) 



= E a ( - K a k 2 t« 



(49) 



Q(x, t\x ) = W(x - x ,t) + W( - x - x , t) , 



(50) 
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where W(x, t) denotes the solution of the FDE (34), for natural boundary conditions. A typical 
example is portrayed in Fig. 8, alongside with its Brownian counterpart in Fig. 9. Similarly, for an 
absorbing barrier, the problem is defined via the Dirichlet condition Q(x , t) — 0, and the half space 
solution for the sharp initial condition W (x) = 8(x — x ) takes on the form 

Q(x, t\x ) = W(x - x , t) - W( - x - x , t) . (51) 

For subdiffusion in a box of extension ( — a, a), the propagator W(x, t) also suffices to determine 
the boundary value problem of two absorbing or two reflecting boundaries which are supposed to 
lie at x — ± a. There, the free solution with the initial value problem W (x) = 5(x) is successively 
folded along the lines through x = ± a, perpendicular to the x axis, i.e., the exact solution is 
constructed with increasing accuracy according to the method of images, to result in the boundary 
value solution [88,205] 

GO 

Q(x, t)= X t w ( x + 4ma > l ) + W(4ma - x + 2a, t)~] , (52) 



m = — oo 



where the minus sign stands for absorbing, the plus sign for reflecting boundaries at x = ± a. We 
note that the solution for the mixed condition of one absorbing and one reflecting boundary is 
obtained via a final folding at the origin of the solution for two absorbing boundaries. Employing 
the relation 13 

a. co / k \ 

£ e~ ikm = e- ikl2 X ( - l) m <5( m + — I , (53) 

m — — co m— co \ / 

we rewrite Eq. (52), after an additional Laplace transform, as 



Q(x,u) = (4a)- 1 X & 



m — — cc 



(54) 



Making use of Eq. (31), the sums can be simplified and evaluated numerically. A typical result is 
shown in Fig. 10, in comparison with the Brownian counterpart. 

An alternative solution procedure, the method of separation of variables is discussed in 
Section 5. 

Being a direct generalisation of Fick's second law, the fractional partial differential equation (34) 
can be solved by standard methods, like the Fourier-Laplace technique, or the method of images. 
The solution W(x, t), Eq. (42), of the FDE (34) has been expressed analytically, in closed form, 



1 3 This relation can be easily proved by the Poisson summation formula [247] 
X f(k)= I dx/(x)e 2 



lu _ L I U -V 

k = — oo m ~ co j oo 

and the integral definition of the delta function [248] 



5( x _ x ') = -L I dfce*"'', 



2n 

as well as some resummations. 
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Q(x,t) 




1 2 3 4 5 

Fig. 9. Brownian image solution Q(x, t) portrayed for the times t = 0.01, 0.1 and 1. The initial condition of starting in 
x = 1 decays much faster in comparison to the subdiffusive case graphed in Fig. 8. 



by the Fox function R\'% which has been studied in detail; compare Appendix B. The pro- 
pagator W(x,t) possesses similar scaling properties like the Gaussian solution (15), i.e., it is 
given by the functional form W(x, t) = (p(t)g(£), where the scaling variable £ is defined through 
£ = x/f lz [249]. 
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Q(x.t) 




Fig. 10. Image solution Q(x, t) for absorbing boundaries at x = +1. Top: The subdiffusive case, a. = 1/2. Bottom: The 
Brownian case, y = 1. Note, again, the cusp shape of the subdiffusive solution. The curves are shown for the times 
t = 0.005, 0.1 and 10 on top and t = 0.05, 0.1 and 10 on the bottom. The broad wings in the top graph are due to the 
stretched Gaussian shape of the free propagator in the subdiffusive regime. 

Verifications of the FDE (34) have so far been discussed for NMR in disordered systems where 
the usually found Gaussian shape of the Fourier transformed propagator gets replaced by an 
inverse power-law pattern [225,250]. These predictions were corroborated by findings from NMR 
experiments in biological tissue [251] which were interpreted along the lines of a fractional 
diffusion equation. 14 A further example is given by FRAP (Fluorescent Recovery After Photo- 
bleaching) experiments [253,254]. 

3.5. Long jumps: Levy flights 



The opposite case of finite characteristic waiting time T and diverging jump length variance 
I 2 can be modelled by a Poissonian waiting time and a Levy distribution for the jump length, i.e., 



X(k) = exp( - o*W) ~ 1 - 0*1*1" 



(55) 



14 It was shown that the results follow an asymptotic power-law in Fourier space which is neither consistent with the 
Gaussian approach nor a Levy flight or walk approach which is of stretched Gaussian shape in Fourier space. 
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for 1 < n < 2, 15 corresponding to the asymptotic behaviour 

M^^A^-^xr 1 -" (56) 

for |x| 5> a. Due to the finiteness of T, this process is of Markovian nature. Substituting the 
asymptotic expansion from Eq. (55) into the relation (25), one obtains 

w ^=vrhw (57) 

from which, upon Fourier and Laplace inversion, the FDE [192,202] 
QW 

— = K»_ a0 D»W(x,t) (58) 

is inferred. The Weyl operator _ ^D^ which in one dimension is equivalent to the Riesz operator V, 
is defined in Section A.2. 16 Here, the generalised diffusion constant is 

K" = (t"/t , (59) 

and carries the dimension [Kf\ = cm' 1 s ~ 1 . The Fourier transform of the propagator can be readily 
computed, obtaining 

W(k, t) = exp( - K"t\k\") , (60) 

which is but the characteristic function of a centred and symmetric Levy distribution, and as such 
used to generate Levy flights [42]. Note that the Fourier space version of Eq. (58) was discussed in 
Ref. [134]. 

In Fig. 11a computer simulation of a Levy flight is shown on the right, in comparison to the 
trajectory of a walk with finite jump length variance I 2 , for the same number of steps. Due to the 
asymptotic property (56) of the jump length pdf, very long jumps may occur with a significantly 
higher probability than for an exponentially decaying pdf like the formerly employed Gaussian 
jump length pdf. The scaling nature of the jump length pdf, as expressed by Eq. (56) leads to the 
clustering nature of the Levy flights, i.e., local motion is occasionally interrupted by long sojourns, 
on all length scales. That is, one finds clusters of local motion within clusters. In fact, the Levy flight 
trajectory can be assigned a fractal dimension d f = fi [38,42,89] and is commonly supposed to be 
an efficient search strategy of living organisms [141, 142,146]. 17 Contrarily, the trajectory drawn 
on the left of Fig. 11, with I 2 < go, fills the two-dimensional space completely, and features no 
distinguishable clusters, as all jumps are of about the same length. 



15 The case < /.i < 1, though very similar, is not discussed here. 

16 Note that we suppress the imaginary unit in Fourier space by adopting the slightly modified definition 
&{-«,D$f(x)} = - \k\" f(k) instead of - m = i"\k\"f(k), following a convention initiated by Compte [192]. 

1 7 Note that we refer to the trajectory of the Levy flight. The spatiotemporal behaviour of these organisms has a finite 
velocity of motion, see below. 



R. Metzler, J. Klafter / Physics Reports 339 (2000) 1-77 



27 




Fig. 11. Comparison of the trajectories of a Brownian or subdiffusive random walk (left) and a Levy walk with index 
ji = 1.5 (right). Whereas both trajectories are statistically self-similar, the Levy walk trajectory possesses a fractal 
dimension, characterising the island structure of clusters of smaller steps, connected by a long step. Both walks are drawn 
for the same number of steps (approx. 7000). 



The solution of the FDE (58) in (x, t) space can again be obtained analytically by making use of 
the Fox functions, the result being [195,217] 



W{x,t) = —H\-\ 



1 

H\x\ 



x 



{KH) 



(UM (1,1/2) 

(1,1), (1,1/2) J 



(61) 



This is a closed-form representation of a Levy stable law, see Appendix C for details. For lim M ^ 2 , 
the classical Gaussian solution is recovered, by standard theorems of the Fox functions. As 
expected, one can infer from Eq. (61) the power-law asymptotics [47,217] 



W(x, t) 



H < 2 



(62) 



typical for Levy distributions. Due to this property, the mean squared displacement diverges: 



(x 2 (t)) -> 00 



(63) 
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which has caused some controversy in that it is different from quantities obtained through cut-offs, 
scaling relations obtained from similarity methods, or rescaling of a fractional moment 



<[x| 5 > oc t' 



din 



(64) 



for < 8 < fi < 2, and loosely called "mean squared displacement" [42,49,52,159,160,227,228]. 
Clearly, Eq. (63) cannot be valid for a particle with a non-diverging mass. For such massive 
particles, a finite velocity of propagation exists, making long instantaneous jumps impossible, see 
below. 

The exact calculation of fractional moments profits from the definition of the Fox functions. 
Thus, writing 



<|x| a > = 2 dxx d W(x,t) , 



(65) 



with < 8, it follows that this integral, due to Eq. (62), only converges for d < \i. In this case, note 
that through Eq. (61) the integral (65) defines the Mellin transformation 



Jt{f{t)\s} 

of the Fox function: 
2 







<M 5 > = 



/"Jo 



2.2 



X 



_(K"t) 



(UM(U/2)" 

(1,1), (1,1/2) . 



2 



H i,i 



2,2 



X 



(K"t) 



Employing the property [237,242] 



.s — 1 um,n 
P>4 



dxx s_1 tf 





(cip, Ap) 


ax 


(b q ,B q )_ 



= a s x(-s) 



of the Fox function, where %(s) is defined in Eq. (B.9), one readily infers 

<\x\ d y = -(K»t?i»Y( -d) = -(K»t) d '» r{ ~ 3 imi + 3) 

W ) ^ t) /( d) t) r{ _ d/2m + g/2) . 



(66) 



(67) 



(68) 



(69) 



Due to the condition < 5 < fi, (\x\ 5 ) is always positive, as can be seen from the T functions in 
Eq. (69). Also, whereas {\x\ 6 } oc t 6 ^ is always sublinear in t, the rescaling of this fractional moment 
leads to the pseudo mean squared displacement [x 2 ] oc t 2 ' 11 , i.e., "superdiffusion". 
Consider two special cases. For 5 -> 0, one proves from Eq. (69) the normalisation: 



lim (\x\ d ) = 1 



(70) 



by help of 1/T(z) ~ z, for z <^ 1. In the Gaussian limit fi = 2, the linear time dependence 



lim <|x|' 5 > = 2K 1 t 



(71) 
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of the mean squared displacement is recovered. Note that numerically obtained fractional 
moments have been used in single molecule spectroscopy for system characterisation 
[138,139]. 

For certain physical systems, for instance the diffusion in energy space encountered in single 
molecule spectroscopy [138,139] or the diffusion on a polymer chain in chemical space [255] the 
divergence (63) of the mean squared displacement does not violate physical principles. However, 
if dealing with massive particles in direct space, physics implies a finite velocity of propagation. 
The latter dilemma can be overcome by replacing Levy flights by Levy walks, a version of 
the CTRW with a spatiotemporal coupling, usually introduced through the delta coupling 
il/(x,t) = iw(t)8(\x\ — vt) [52,151-153], with a time cost penalising long jumps [49,52,151]. This is, 
for example, the case in the chaotic phase diffusion in Josephson junctions [256], in turbulent flow 
[257], or in chaotic Hamiltonian systems [49,258,259]. As we are mainly concerned with subdif- 
fusive systems, Levy walks are not considered here. 

Boundary value problems for Levy flights are more involved than for the subdiffusive case, as 
the long jumps make the definition of a boundary actually quite intricate. Such problems were 
discussed for Levy flights in the half-space by Zumofen and Klafter [260] and for a box by 
Drysdale and Robinson [261]. 



3. 6. The competition between long rests and long jumps 

As noted, the problem of the diverging mean squared displacement encountered in the discussion 
of Levy flights is often circumvented by the consideration of (x-t) scaling relations, or measuring 
the width of the pdf W(x, t) rather than its variance. An alternative method was applied in 
Ref. [217], making use of the definition 



<X 2 (f)>z 



dxx 2 W(x,t) ~ t 2l » (72) 



according to which the walker is considered in an imaginary, growing box, leading to the ~ t 2,,i 
behaviour. The latter was verified by numerical simulations [217]. Note that the cut-offs of 
the integral in Eq. (72) are time dependent. The imaginary box spans the spatial interval 
A(t) = (L 1 — L 2 )t xl!1 which grows in the course of time. It gives a measure, that a finite portion of 
the probability is gathered within the given interval A(t). 

By use of such relations, one can consider a random walk characterised through broad pdfs for 
both waiting time and jump length, thus leading to infinite T and I 2 , and the FDE 

dW 

— = D}-'KWW(x,t) (73) 

with K£ = a^jx 1 . In this case, the appropriately defined quantity <x 2 (r)> L , which we call the pseudo 
or imaginary mean squared displacement, reveals the temporal form 



<X 2 (f)>L ~ t 2 ^ 



(74) 
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Fig. 12. Phase diagram for the FDE (73). The different "phases" characterise the four domains which can be distin- 
guished according to diverging or finite characteristic waiting time T and jump length variance Z 2 . If the parameters 
a and p. become larger than 1 and two, respectively, the corresponding CTRW processes lock onto processes dominated 
by finite T or Z 2 . For Z 2 -> oo , one observes either the traditional Markovian Levy flights for a > 1, or the 
non-Markovian version (s) with single modes decaying in a Mittag-Leffler pattern. 



Consequently, the phase diagram displayed in Fig. 12 can be drawn. Note that in expression (74), 
the exponent a enters in a "proper" way, i.e., it is already present for I 2 < go, whereas the Levy 
index \i manifests the artifice of introducing the pseudo mean squared displacement. Conversely, 
the competition of laminar motion events ("flights") and localisation (waiting) events in the Levy 
walk picture is given through the relation [262] 



whereby the time spans spent in laminar motion events is governed by the pdf ij/(t) ~ at' 11 ' \ and 
the waiting times are drawn from the pdf \j/(t) ~ bt'*' 1 . Eq. (75) for Levy walks clearly differs from 
the multiplicative nature by which the exponents a and \x enter into Eq. (74). 

3. 7. Whafs the course, helmsman! 

Let us take stock at this point and address two important issues concerning the fractional 
diffusion concept. 

3.7.1. The long time limit and its consequence for the fractional diffusion equation 

In the derivation of both the subdiffusive FDE (34) and the Levy flight FDE (58) from the CTRW 
scheme, the diffusion limit (k, u) -> (0, 0) was drawn. Equivalently, the diffusion limit corresponds 
to choosing (a, t) -> (0, 0) with K a = <t 2 /t sc = const or K 11 = a^/t = const which matches the limit 




(75) 
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K = lim Ax ^o,Ai^o(Ax) 2 /(2At) drawn in the Brownian case. In this sense, the FDE is valid in the 
diffusion limit t > t. In this diffusion limit, the equivalence to the CTRW model holds true for all 
moments. Indeed, from the FDE (34) one obtains the general expression 

<x 2 "(t)> = (2n)\ " . (76) 
1(1 + not) 

From the CTRW theory with untruncated wave number expansion, one finds higher-order 
corrections of the form 

<x 2 "(f)> = (2n)! ' (1 + OdKjT 1 )) (77) 
1(1 + na) 

which can be neglected in the drawn diffusion limit. 

Differences with the CTRW model occur in details of the pdf W(x, t) where higher orders in the 
k expansion come into play [239], and in higher dimensions as discussed in Ref. [215]. 



3.7.2. Fractional diffusion equations and limit theorems 

The Gaussian solution obtained for Brownian motion is a consequence of the central limit 
theorem. How should the results for the subdiffusive FDE (34) and its Levy flight counterpart (58) 
be judged from the limit theorem point of view? 

The subdiffusive process combines the long tailed waiting time process with a jump length 
distribution that possesses a finite characteristic variance. The resulting pdf W(x, t) which is the 
solution of the FDE (34) is not Gaussian. However, in the sense of Bouchaud and Georges [42], 
there exists a central limit theorem controlling this process as W(x, t) fulfils a scaling relation 
W(x, t) = ^(t)f(x/^(t)) as can be seen from Eq. (42). Moreover, the temporal distribution of this 
subdiffusive process is a one-sided Levy distribution characterised through its Laplace transform, 
the quantity W(x, u) as given in Eq. (44). The existence of such a central limit theorem for this type 
of random walk process is further corroborated by the properties of the Fox function in the 
solution (42) which guarantee a smooth transition to the classical Gaussian in the limit a -> 1. 

Levy flights are Markovian processes and they are governed by a pdf which is Levy stable. 
Therefore, they are a direct result of the Levy-Gnedenko generalised central limit theorem. 

Thus, the FDE concept is a valid extension of the standard diffusion equation in the same sense 
as CTRW generalises Brownian motion. 



4. Fractional diffusion-advection equations 

Diffusion with an additional velocity field v and diffusion under the influence of a constant 
external force field are, in the Brownian case, both modelled by the diffusion-advection equation 
(DAE) 



dW ^dW _ 8 2 
Qt 8x 1 dx 



v— = K 1 ^W(x,t) . (78) 
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In the case of anomalous diffusion this is no longer true, i.e., the fractional generalisation may 
be different for the advection case and the transport in an external force field. We start 
with considering the external velocity field, the external force field is discussed in the following 
section. 

In the stationary state, i.e., without inertial terms, described by the DAE (78), the problem is 
Galilei invariant: it is invariant under a transformation x -> x — vt. Requiring this property to carry 
over to the anomalous case, a straightforward extension of the CTRW scheme leads to a fractional 
diffusion-advection equation (FDAE). A Galilei variant model is discussed subsequently. Alterna- 
tive formulations for advected Levy flights are presented at the end of this section. 

A more detailed summary of FDAEs including dispersive sedimentation processes and the 
partial sticking mechanism, as well as suggestions for the measurement, are given in Ref. [207]. 



4.1. The Galilei invariant fractional diffusion-advection equation 

In the moving frame of the test particle which is dragged along the homogeneous velocity field v, 
the random walk is governed by the usual jump pdf \j/(x, t). The corresponding jump pdf <p(x, t) in 
the laboratory frame is consequently obtained via the transformation 

4>(x, t) = \j/{x - vt, t) . (79) 

This carries over, Fourier-Laplace transformed, to the functional relation 

4>{k, u) = i]/(k, u + ivk) (80) 

between <p(k, u) and \j/{k, u). For the case of infinite characteristic waiting time T, i.e., a broad 
waiting time pdf, and finite jump length variance I 2 , the Fourier- Laplace form of the propagator, 



1 

u + ivk + KJi 2 u l 



W ( k > ") = „ ,, : ; , (81) 



can be directly deduced from the CTRW solution (25) for the diffusion limit k -» and u -> 0. 
Proceeding along the same steps as outlined in Section 3, one is led to the FDAE [202] 

— + ,— - M -K. w WM, (82, 

by virtue of Eq. (32). Due to the required Galilei invariance, the solution for the propagator in (x, t) 
space is given by the Galilei-shifted solution of the FDE (34), i.e., 

W(x,t)= W v=0 (x-vt,t) (83) 

where W v = (x, t) denotes the free propagator according to Eqs. (42) and (43), for the sharp initial 
value W (x) = d(x). 
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The moments of the FDAE (82) are readily calculated, and one obtains 

<x(t)> = vt , (84) 

2 K 
1 (1 + a) 

<(Ax(t)) 2 >= — -f-^t". (86) 
1 (1 + a) 

Thus, the mean squared displacement <(Ax(t)) 2 > contains solely the molecular contribution, i.e., the 
relative mixing in the moving frame, whereas the first moment <x(£)> accounts for the simple drag 
along the velocity field v, as is to be expected from the ordinary drift term v 8 W/dx occurring in the 
FDAE (82). The assumed Galilei invariance thus carries over to the moments, as it should. This 
behavior is also found in the temporal evolution of the pdf which is depicted in Figs. 13 and 14. 

A possible realisation of the Galilei invariant subdiffusion is the motion of a particle in a flow 
field where the flowing substance itself causes the occurrence of subdiffusion, like in the case of 
a polymer solution, and a small bead immersed in it. 

4.2. The Galilei variant fractional diffusion-advection equation 

Instead of assuming the Galilei invariance expressed in Eqs. (79) and (80) leading to the 
propagator (83), Compte [205], and Compte et al. [206], for an explicitly position-dependent 
velocity field v(x), assume the relation 

(j)(x, t; x ) = <A(x - T a v(x ), t) (87) 




Fig. 13. Galilei invariant subdiffusive model. The propagator is shown for the dimensionless times t = 0.02, 0.2 and 2. 
The propagator is symmetric with respect to its maximum which is translated with velocity v = 1. The cusps marking the 
initial condition are distinct in comparison to the Brownian result shown in Fig. 14. 
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Fig. 14. Galilei invariant Brownian model. The propagator is shown for the dimensionless times t = 0.02, 0.2 and 2. The 
propagator is symmetric with respect to its maximum which is translated with velocity v = 1. The fast decay of the initial 
condition is mirrored in the smooth peak of the distribution, compare to the subdiffusive result shown in Fig. 13. 



between the jump pdf (/> and the jump pdf \jj(x, t) of the free diffusion process. Thereby, </> depends 
on both the jump length x and the starting point x . Additionally, a microscopic advection time 
T a is introduced. On that basis, the framework developed in Ref. [205] is applied in Ref. [206] to 
Taylor flows. 

The resulting FDAE [204-206] 



QW 
~d7 



1 -a 



- A a -v(x) 



K 



a 8x 2 



W(x,t) 



(88) 



has the same structure as the fractional Fokker-Planck equation. For a homogeneous velocity 
field, the resulting FDAE 



QW 

~0T 



= oA 1 



„ 8 v 82 ' 
- A —v + K a 

ox ox 



W(x, t) 



(89) 



has a constant drift coefficient. It can be proved that the fractional solution does not fulfil 
a generalised Galilei invariance of the form W(x — v*f, t). However, according to Section 5.4 the 
fractional solution W x (x, t) can be expressed in terms of the Brownian solution W^x, t), the Galilei 
shifted Gaussian 



W 1 (x,t) = 



1 



AnKt 



exp 



(x - vtf 
4Kt 



(90) 



for numerical purposes [207]. This way, Fig. 15 was obtained which is to be compared with Figs. 13 
and 14 for the Galilei invariant case. Note the growing skewness of the solution. This persistence of 
the initial condition has already been celebrated by Scher and coworkers. 
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Fig. 15. Galilei variant subdiffusive model. The propagator is shown for the dimensionless times t = 0.02, 0.2 and 2. The 
propagator is asymmetric with respect to its maximum which stays fixed at the origin. The plume stretches more and 
more into the direction of the velocity. 



The difference to the Galilei invariant model is also mirrored in the moments. From Eq. (89) one 
deduces 



<x(t)> = 



A a vf 

r(i + a) 



2A 2 v 2 t 2a 



<x 2 (t)> = 



+ 



2K a f 



r(l + 2a) T(l + a) 



(91) 



(92) 



where now the first moment increases sublinearly in time. This model describes physical systems 
where trapping occurs, compare Section 6, i.e., the particle gets repeatedly immobilised in the 
environment for a trapping time drawn from the waiting time pdf w(t), before it gets dragged along 
the velocity stream again. Experimental realisations might be found in porous systems like the one 
displayed in Fig. 16 where the particle gets trapped in still regions off the velocity backbone, or the 
multiple trapping systems addressed in Section 6, or possibly in gel electrophoresis [263]. 

4.3. Alternative approaches for Levy flights 



For Levy flights in an external velocity field v, i.e., for finite characteristic waiting time T but 
infinite jump length variance I 2 , the FDAE [202] 



dW dW 

+ v— = XW(x,t) 



dt 



Qx 



(93) 



is recovered [202], which describes a Markovian process with diverging mean squared displace- 
ment, see the discussion in Section 3.5. The drift term vdW/dx in Eq. (93) is the same as in the 
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Fig. 16. Swiss cheese percolation network, as used by Klemm et al. [106]. The image clearly shows the relatively sharp 
distinction between regions with high flow velocities (backbones), and still domains. A particle trapped in one of the 
side-pores off the backbones stays approximately stationary, until it rejoins the main stream backbones. Courtesy 
A. Klemm and R. Kimmich, Ulm University. 



standard DAE (78), and thus the result is simply given in accordance to the subdiffusive case, 
Eq. (83), by W v (x, t) = W p = (x — vt,t) where here W v = (x, t) is the Levy stable solution (61). The 
latter is symmetric in x. 

There may occur situations where the resulting propagator is asymmetric. This case is modelled 
by the alternative approach in Refs. [208,209] where a skewed Levy distribution is assumed instead 
of the symmetric one leading to the Riesz/Weyl operator V introduced in the preceding parts. 
Consequently, the resulting FDAE possesses an additional parameter accounting for the asym- 
metry, in comparison to the FDAE (93). 



5. The fractional Fokker-Planck equation: anomalous diffusion in an external force field 

Many physical transport problems take place under the influence of an external force field: 
a constant electrical bias field exerting a force on charge carriers, a periodic potential encountered 
in certain problems in solid state physics or in the modelling of molecular motors, a bistable 
potential in reaction dynamics or molecular switching processes, or a harmonic potential describ- 
ing a bound particle. In this section, a framework for the treatment of anomalous diffusion 
problems under the influence of an external force field is developed. A physical model based on the 
Langevin equation with Gaussian white noise for multiple trapping systems is introduced in the 
next section. 
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5.1. The Fokker-Planck equation 



Normal diffusion in an external force field is often modelled in terms of the Fokker-Planck 
equation (FPE) [24,29,35-37,78-80,264-266] 

dW 



8 V'(x) 

+ ^l^"2 



8x mr\x dx z 



W(x, t) (94) 



where m is the mass of the diffusing test particle, n 1 denotes the friction constant characterising the 
interaction between the test particle and its embedding, and the force is related to the external 
potential through F(x) = — dV(x)/dx. Let us first list some basic properties of the FPE (94) to 
which the fractional counterpart developed below can be compared. 

(i) In the force-free limit, the FPE (94) reduces to Fick's second law and thus the time evolution of 
the mean squared displacement follows the linear form (4). 

(ii) Single modes of the FPE (94) relax exponentially in time, 

r n (0 = exp(-A n>1 0, (95) 

where A n>1 is the eigenvalue of the Fokker-Planck operator L FP defined below. 

(iii) The stationary solution 

W st (x) = lim W(x, t) (96) 

t->QO 

is given by the Gibbs-Boltzmann distribution 
W st (x) = N exp( - pV(x)) (97) 

where JV is a normalisation constant and /? = (/c B T) _1 denotes the Boltzmann factor. 

(iv) The FPE (94) further fulfils the Einstein-Stokes-Smoluchowski relation [38] 

K^kvT/mr,, (98) 

which is closely connected with the fluctuation-dissipation theorem. 

(v) Finally, the second Einstein relation 

<«„> ,99) 
2 /c B T 

is recovered from Eq. (94) which connects the first moment in presence of the constant force F, 
(x(t)) F , with the second moment in absence of this force, 

<x 2 (t)> =2K,t . (100) 
The latter relationship, Eq. (99), is a consequence of linear response. 



5.2. The fractional Fokker-Planck equation 



The FPE (94) is well studied for a variety of potential types, and the respective results have found 
wide application. For the description of anomalous transport in the presence of an external field, 
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we introduce a fractional extension of the FPE, namely the fractional Fokker-Planck equation 
(FFPE) [212-215] 



m - D 1 



8 V'(x) | ^ 8 2 
8x mn a *9x 2 



W(x,t) . (101) 



The FP-operator 

9 V'(x) | 6 2 
8x m?; a a 9x 



occurring in the FFPE (101) contains the generalised diffusion constant and the generalised 
friction constant ^ of dimension [f/ a ] = s*~ 2 . For a -> 1, the standard FPE (94) is recovered, and 
for V(x) = const., i.e. in the force-free limit, the FDE (34) emerges. 
We will show that this equation fulfils the following requirements: 

(i) The FFPE (101) describes subdiffusion in accordance to the mean squared displacement 

2K 

^-mt? (103) 

in the force-free limit. This is obvious as for V(x) = const., Eq. (101) reduces to the standard 
FDE, Eq. (34). 

(ii) The relaxation of single modes is governed by a Mittag-Leffier pattern 

T n (t) = E a ( - X nA f) . (104) 

(iii) The stationary solution is given by the Gibbs-Boltzmann distribution (97). 

(iv) A generalisation of the Einstein-Stokes-Smoluchowski relation (98) connects the generalised 
friction and diffusion coefficients. 

(v) The second Einstein relation (99) can be shown to hold true for Eq. (101). 

The FFPE (101) will be derived explicitly in the following subsection, and independently in the 
next section as the high friction limit of a fractional phase space model. 

In order to derive the stationary solution W st (x) of the FFPE (101), note that the right-hand side 
of the Eq. (101) can be rewritten as 

in terms of the probability current [36] 

SM = {-^- K i) WM - (106) 

If a stationary state is reached, S(x, t) must be a constant. Thus, if S st (x ) = at any point x , it 
vanishes everywhere, and the stationary solution of the FFPE satisfies [36,212] 

^-WJx) + K x ^-WJx) = (107) 



R. Metzler, J. Klafter / Physics Reports 339 (2000) 1-77 



39 



from which the exponential result 

W st (x) = Ncxp( -"^-) (108) 

can be inferred. Requiring, in analogy to the standard case, that W st is given by the Boltzmann 
distribution (97), the generalised Einstein-Stokes-Smoluchowski relation [42,203,212-215,267] 

K x =k B T/mn x (109) 

is readily recovered. Thus, the FFPE (101) obeys some generalised fluctuation-dissipation the- 
orem. The generalised Einstein-Stokes-Smoluchowski relation (109) has recently been supported 
by findings of Amblard and coworkers [120,121]. In order to check whether the FFPE (101) 
satisfies the second Einstein relation (99), the first moment in presence of the constant force F is 
calculated under linear response conditions, obtaining the expression 

<x(t)> F = — Jrr—-/- ( n °) 

mt] a r(l + a) 

A comparison to the force-free result of the mean squared displacement, Eq. (103), leads to the 
second Einstein relation (99). The validity of this relation in the case of subdiffusion has found 
support experimentally in the work by Schiff et al. [103]. 

At this point it is worth digressing a bit and noting that the condition dW/dt = 0, i.e., the usual 
stationary condition, allows for a unique solution for the fractional equation (101), namely 
the Gibbs-Boltzmann form (97). On the one hand, dW/dt = implies W = const(x). On the 
other hand, D}~ a W(x, t) = is fulfilled by W(x, t) = const(x) corresponding to the reasoning 
in the above derivation via the vanishing probability current, S sl (x) = 0, as well as by 
W(x, t) = const(x)t~ a , according to relation (37). This contradicts, however, the requirement 
W — const(x) imposed by the stationarity condition. Thus, the exponential form derived above is 
unique. This fact is mirrored in the stationary solution obtained by the method of separation 
of variables discussed below, which leads to a time-independent solution. Note further that a quasi- 
stationary form similar to xf 1 was found for a fractional phase space equation by Hilfer in his 
thermodynamical derivation of fractional dynamics [235,268,269]. 

5.3. Separation of variables and the fractional Ornstein-Uhlenbeck process 

5.3.1. The separation of variables 

In order to determine the pattern according to which the stationary solution is approached, the 
separation ansatz 

W n (x,t) = T„(t)(p„(x) (111) 

for a given mode n is introduced into the FFPE (101), leading to the decoupled set of eigen- 
equations 

-A n>a oD}-«T n (t) (112) 



dr„(t)_ 



dt 



L FP (p n (x) = 



(113) 
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featuring the fractional eigenvalues A„ >a . The temporal eigenfunction T n (t) being governed by the 
fractional relaxation equation (112), it is described by the Mittag-Leffler pattern 

T M = E.(-X„,.f)=Y^l (114) 

with the choice T„(0) = 1. For a = 1, the standard exponential form T„(t) = exp( — l„^t) follows, 
whereas for < a < 1, the initial stretched exponential behaviour 

TM ~™i~Mi) ,115 » 

turns over to the power-law long-time behaviour 

T -«>~nd*u?- ,116) 

This interpolation from stretched exponential to inverse power-law behaviour has been reported 
from the rheology of polymeric systems [20,21], and from protein rebinding [22]. Let us briefly 
examine the importance of the Mittag-Leffler function. In Refs. [180,270] it is shown, that the 
Mittag-Leffler function is the exact relaxation function for an underlying fractal time random walk 
process, and that this function directly leads to the Cole-Cole behaviour [270,271] for the complex 
susceptibility which is broadly used to describe experimental results. Furthermore, the Mit- 
tag-Leffler function can be decomposed into single Debye processes, the relaxation time distribu- 
tion of which is given by a modified, completely asymmetric Levy distribution [270]. This last 
observation is related to the formulation of Mittag-Leffler relaxation described in Refs. [21,180]. 
In Ref. [245], the significance of the Mittag-Leffler function was shown, where its Laplace 
transform was obtained as a general result for a collision model in the Rayleigh limit. 
The full solution of the FFPE (101) is given by the sum over all eigensolutions, i.e., by 



W(x, t\x',0) = e^'V 2 -*w/ 2 X U*)Ux')E y ( - X n ,X) (117) 

n 

for an initial distribution concentrated in x'. Here, the eigenfunctions \j/ n {x) = Q^^cp^x) of the 
Hermitian operator L are defined in terms of the <p n (x), the eigenfunctions of the FP-operator L FP , 
via the scaled potential <t>(x) = F(x)/[/c B T], L and L FP sharing the same eigenvalues A„ >a [36]. On 
arranging these eigenvalues in increasing order, i.e., < X 0<a < A 1>a < A 2 ,« < the first eigen- 
value is zero iff there exists a stationary solution, which is positive definite and defined through 
W si (x) = lim,.^ W(x, t). This stationary solution is independent of the fractional index a and 
coincides with the required Boltzmann solution (97). However, the temporal relaxation of a single 
mode n towards equilibrium is highly non-exponential. 

The FFPE (101) obeys, in Laplace space, the functional relation [212] 



H^x, li) — ——u x 1 W± ( x, ——u x j 



(118) 
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for the same initial condition W (x) = S(x — x'). The subscripts refer to the solutions of the FPE 
(94) for a = 1 and the FFPE (101) for arbitrary < a < 1, respectively. Thus, in the Laplace 
domain, the subdiffusive system passes through the same states as its normal counterpart, for 
a rescaled Laplace variable. In Section 5.4 we show that the solution W x of the FFPE (101) can be 
expressed in terms of Wx through an integral relation from which the non-negativity of W a can be 
proved, as the Pawula theorem guarantees that the solution W l is a proper pdf. For the force-free 
FDE (34) describing subdiffusion, Schneider and Wyss proved directly that the propagator W(x, t), 
Eq. (42), is a proper pdf. 



5.3.2. The fractional harmonically bound particle 

It is worthwhile considering the example of a subdiffusive, harmonically bound particle, i.e. the 
subdiffusive motion in the potential V(x) = jmco 2 x 2 which exerts a restoring force on the test 
particle. We find the following solution [212] 



employing reduced coordinates t = tjx and x = Xy/ma> 2 /[k B T~\, as well as x~ a = co 2 /ri a . H n denotes 
the Hermite polynomials [236], and the eigenvalues here are A„ >a = nco 2 /ri a . This result is plotted in 
Figs. 17-19 in the course of time, for an asymmetric initial condition, and compared to the 
Brownian result. The stationary solution of this process is found to be [236] 



i.e., the Gibbs-Boltzmann distribution, as it should. 

For a given potential V(x), one can extract the moments of order n, <x"(t)>, directly from 
the FFPE (101), by integration over J*^ dxx". This leads to a fractional differential equation 
for the respective moment. For instance, for <x(f)> one obtains the relation (d/df)<x(f)> 
= oA 1 ^ft^A/aK^WX i- e -> the fractional relaxation equation, solved by the Mittag-Leffler pattern 

<x(0> = x E a ( - (t/xf) . (121) 

Eq. (121) describes the relaxation of the mean of an off-centre initial distribution, sliding into the 
symmetric final state on the bottom of the potential valley, which is characterised by the symmetric 
limit lim,^ Q0 <x(t)> = 0. The temporal evolution of the second moment, 

<x 2 (f)> = xf h + [_x 2 - xf„-]E a ( - 2(t/xY) (122) 

also follows the Mittag-Leffler pattern, however characterised by the relaxation time scale t/2 1/=c . 
Hereby, the thermal equilibrium is defined as x 2 h = k B Tl\ma) 2 ~\ reached for t -> go . The second 
moment (122) is graphed for the subdiffusive case in Fig. 20, in comparison to the Brownian 
counterpart. 

Note that this process discussed in position space is equivalent to the fractional version of the 
Ornstein-Uhlenbeck process. The Laplace transformed versions of Eqs. (121) and (122) were found 
in Ref. [245] to describe the relaxation of a heavy test particle immersed in a bath of light particles, 
in a generalised Rayleigh limit. 
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W(x,t) 




Fig. 17. Pdf W(x, £), Eq. (119), of the fractional Ornstein-Uhlenbeck process, for the anomalous diffusion exponent 
a = 1/2. The initial value is chosen to be W (x) = 5(x — 1). The maximum clearly slides towards the origin, acquiring an 
inversion symmetric shape. The curves are drawn for the times t = 0.02, 0.2, and 20, employing the integral relation with 
the Brownian solution. Note the distinct cusps around the initial position. Compare Fig. 18. 



W 




Fig. 18. Fractional Ornstein-Uhlenbeck process. Comparison of the numerical behaviour of the summation representa- 
tion (dashed) with 151 summation terms and the integral representation (A transform). The latter is obtained by 
Mathematica employing the numerical integration command Nlntegrate. The cusp which is a typical feature for 
subdiffusive processes is much more pronounced in the curve obtained through the integral transformation. The 
computation time for the latter is even shorter than for the calculation of the truncated sum so that this representation is 
preferable for numerical purposes. 
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-4-2 2 4 

Fig. 19. Fractional Ornstein-Uhlenbeck process. Comparison of subdiffusive pdf at t = 0.4 and 20 with Brownian 
solution at t = 0.36 (dashed) and the stationary solution (dashed). The cusp emphasising the slowly devaying initial 
condition at x = 1 in the subdiffusive case is distinct. Also note the pronounced asymmetry of the fractional solution. 



5.4. The connection between the fractional solution and its Brownian counterpart 

The solution of the FFPE, W(x, t), follows the scaling relation given in Eq. (118). As shown by 
Barkai and Silbey [272], Eq. (118) can be rewritten in the form 

W x {x,t)=\ dsA^^W^s) (123) 

which corresponds to a modified Laplace transformation from t to {r\^r\ x )u^-. The kernel A(s, t) is 
defined in terms of the inverse Laplace transformation 



the result being the modified one-sided Levy distribution L a 



+ 18 



1 - +/ * 



A(s,t) = ), * = whi. (125) 

Consequently, the transformation (123) guarantees the existence and positivity of W x (x,t) if only 
the Brownian counterpart, Wi(x, t), is a proper pdf. 



18 L a + (f) is defined on the positive real axis and shows the asymptotic power-law behaviour L a + ~ t 1 a , compare 
Appendix C. 
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Fig. 20. Mean squared displacement for the fractional (a = 1/2, full line) and normal (dashed) Ornstein-Uhlenbeck 
process. The normal process shows the typical proportionality to t for small times, and approaches the saturation value 
much faster than its subdiffusive analogue which starts off with the t 1/2 behaviour and approaches the thermal 
equilibrium value by a power-law, compare Eq. (122). 

In the normalised version, the kernel A is given through A(s, u) — u a ~ 1 e~ su ". Rewriting it in terms 
of a Fox function, one obtains the exact representation 



s l/« 


(1,1) " 


t 


(1, !/«)_ 



(126) 



in terms of the Fox function if};? whose series representation reads 



From the properties of the Fox function, or by using Mathematica, one can derive the following 
results for some special cases: 




(127) 



a = 1/2: 




(128) 



a = 1/3: 




(129) 
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a = 2/3: 




1 1 

^|_r(T/3) 




m 2 



4s 3 



t 2/3 r( - 1/3) 



s 




(130) 



With such expressions at hand, one has a very convenient method for plotting the fractional 
solutions via relation (123). This method shows a good performance when used in Mathematica 
with the NIntegrate command. Via this method, Figs. 15 and 17-19 have been obtained. 

Note that the connection between 147, and W 1 is related to the quantity x s (t) that exactly s events 
occurred in time t known from a random walk characterised by the waiting time distribution 
w(t), Xs( u ) = E W ( M )] S [1 — w(w)]/w(w) [277]. For the special form of the waiting time distribution, 
w(m) = e~"", one finds in the long time limit x s (u) = u*~ 1 e~ su * = A(s, u). 

5.5. The fractional analogue of Kramers escape theory from a potential well 

As we have seen, it is an important feature of the fractional kinetic equations that its solution 
W x can be expressed in terms of its Brownian counterpart W-^. In fact, all kinetic processes 
associated with such a fractional equation are affected by scaling relations such as Eq. (118). 

Let us recall that in the standard Kramers problem the escape of a scalar test particle subject to 
a Gaussian white noise over a potential barrier is considered in the limit of low diffusivity, AV/K 
where AV is the barrier height and K the diffusion constant [36]. The temporal decay of the 
probability to still find the particle within the potential well is given by an exponential function 



with AV = F(x max ) — F(x min ); compare Fig. 21. In Eq. (132), the exponential function contains the 
Boltzmann factor ft = (k B T) ~ 1 so that the inverse Kramers rate follows an Arrhenius activation 
fK 1 oce c/r [274]. 

Let us now derive the fractional counterpart to the exponential decay pattern from Eq. (131). 
Application of relation (118) to the Laplace transform p(u) = (r k + u)' 1 of the survival probability, 
Eq. (131), produces 



p(t) = e 



(131) 



where the Kramers rate in the overdamped limit is defined through [36] 




(132) 



pM = 



(133) 



with the generalised, fractional Kramers rate 




1 



V F "(x mi „)IHw)|exp( - P AV) . 



(134) 



2nmrj 



a 
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V(x) 




Fig. 21. Potential well in the Kramers rate model. Initially the particle is assumed to be caught in the potential hole. The 
x axis corresponds to a reaction coordinate. 



Note that consequently the Arrhenius form of the temperature activation is preserved. Via Laplace 
inversion of Eq. (133) one finds the fractional survival probability 

pM = E x ( - r£Y) (135) 

in terms of the Mittag-Leffler function E a . 

Often the notion of time-dependent rate coefficients is preferred, i.e., the survival probability is 
defined as p(t) — exp( — k(t)t) in terms of the rate coefficient k(t). For the fractional Kramers model 
we therefore find k(t) = \lnE x ( — r$t*)\/t which leads to the two limiting cases 

< 136 > 

and 

k(t) ~ yln(t[r^r(l - a)] 1/a ), t (r^) 1/a . (137) 

A detailed investigation of this topic is found in Ref. [273], including the discussion of a possible 
application to ligand rebinding in proteins, compare Refs. [200,275,276]. 

5. 6. The derivation of the fractional Fokker-Planck equation 

For the derivation of the subdiffusive FFPE (101), and a further generalisation accounting 
for a broad jump distance statistics introduced below, we again draw from the random walk 
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formulation and its duality with (generalised) diffusion equations. In order to derive the standard 
FPE, the random walk model employed in Section 3 has to be modified as to account for the 
broken spatial homogeneity, in comparison to the force-free diffusion under natural boundary 
conditions. A non-linear external force field F(x) acting upon the system, the local probabilities to 
jump right or left, A(j) and B(j), explicitly depend on the position j. The corresponding discrete 
master equation therefore becomes 

Wj(t + At) = Aj. t Wj-.it) + B j+ ! W j+1 (t) , (138) 

compare to Eq. (10). Taylor expansions in time, analogous to those already introduced, and in 
space, 

, „ 7 „ x a QA(x)W(x,t) (Ax) 2 Q 2 A(x)W(x, t) _ /r4n ,, /1om 
Aj^Wj-M = A(x)W(x,t) - Ax 8x + ~y Qx 2 + 0([Ax] 3 ) (139) 

lead to the FPE (94), with the appropriate limits 
V'(x\ Ax 

— = lim —lB(x)-A(x)-], (140) 

Kl = lim (I4l) 

For taking these limits, we impose the normalisation A(x) + B(x) = l and note that the in- 
homogeneity in jumping left or right, A(x) — B(x), becomes small for Ax -> 0, according to 
a Boltzmann distribution for a system close to thermal equilibrium, so that 

B(x) - A(x) ~ + 0([Ax] 2 ) . (142) 

The master equation (138) involves local steps in time and position, thus enabling the Taylor 
expansions. Broad waiting time pdfs or jump length pdfs do not allow for such an expansion, 
as they are connected with long-range steps, and thus Ax or At cannot be considered small 
parameters. 

An extension of above scheme for subdiffusion in an external force field has recently been 
presented by Barkai et al. [215], and is discussed below. Here, we prefer an alternative derivation 
which also accounts for a broad jump length statistics [213,214]. Note that for non-local jumps, 
already the simple master equation (138) has to be modified to the form 

Wj(t + At)= t A ln Wj- n (t)+ I B ln W j+n (t) (143) 

n=l n=l 

allowing for jumps from any site j ± n to site;', involving the normalisation 

00 

X (A hn + B hn ) = 1 , (144) 

n=l 
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i.e., a step to site j can come from any site j ± n. In order to introduce the continuum limit of this 
random walk model, the transfer kernel [213,214] 



A(x,x') = X{x - x')[A(x')0(x - x') + B(x')0{x' - x)] 



(145) 



replaces the homogeneous and isotropic jump length pdf X(x — x') = X(\x — x'|) in the CTRW 
approach. Thus, the function A explicitly depends on both the departure site x' and the arrival site 
x. A obeys the normalisation condition 



ddA(x'\\8) = 1 



(146) 



where /l(x'||x — x') = A(x, x'). Including an arbitrary waiting time pdf w(t), it was shown that the 
underlying inhomogeneous CTRW is governed by the generalised master equation [213] 



QW 

with the kernel 



dx' dt' K(x, x'; t - t')W(x', t') 



K(x, x'; u) = uw(u) 



A(x, x') — 5(x) 



1 — w(u) 

given in Laplace space; or, equivalently, through the equation 



W(x, t) 



dx' 



dt'w(t - t')A(x,x')W(x',t') + W(t)W (x) . 



(147) 



(148) 



(149) 



Transforming to Fourier-Laplace space, the FFPE (101) can be derived for a broad waiting time 
pdf with a diverging characteristic waiting time T, Eq. (20), in combination with a finite jump 
length variance I 2 , Eq. (21). Taking also non-local jump statistics into account, i.e., assuming 
a jump length pdf with infinite I 2 , one recovers the FFPE 



QW 



oA 1 



8 V'(x) 



W(x, t) 



dx mt] a 

whereby the drift and diffusion coefficients are given by [213,214] 

mrj a /it 



(150) 



(151) 



Kit 



(152) 



The FFPE (150) thus describes the competition between subdiffusion and Levy flights, as it was 
already encountered for the FDE (73). 

An alternative derivation for the FFPE (101), from an asymmetric random walk scheme is given 
by Barkai et al. as follows [215]. 



R. Metzler, J. Klafter / Physics Reports 339 (2000) 1-77 



49 



We consider an unbounded random walk on a one-dimensional lattice with a lattice spacing a. 
Individual lattice sites are denoted by { . . . , — 1, 0, 1, . . . , n, . . . }. At time t = the particle be located 
at site n = 0. Once the particle has arrived at site n it is trapped there for some random time. These 
waiting times are given by {tJ and = 1,2,.... The {tJ are independent random variables, 
identically distributed according to a pdf w(t). It is assumed that w(t) is independent of the location 
of the particle n (i.e., it is independent of the external field). It is further assumed that the particle 
executes only nearest neighbour jumps. The probability of hopping from site n to n + 1 is A(n), and 
from site n to site n — 1 it is B(n), the normalisation condition being A(n) + B(n) = 1. A(n) and B(n) 
are time independent. 

The probability that the random walker has jumped i times in the interval (0, t), Qi(t), is given in 
Laplace space through 

Qi(u) = l -^^w(uY , (153) 
u 

using the Laplace transform W(u) = (1 — w(u))/u of the sticking probability (24). If W„(t) denotes the 
probability of finding the particle at site n at time t, and Pi (n) be the probability to be on site n after 
step i, then 

W„(t) = £ Pi(n)Qi(t) . (154) 



i = 



Using Eq. (153), 

W n {u) = ^ X Pi (n)w(uy . (155) 

u i = 

The evolution of pi(n) is determined by the discrete time and space equation 

p i+1 (n) = A(n - \) Pi (n - 1) + B(n + l) Pi (n + 1) . (156) 

In Eq. (156) we have used the assumption that the directional jump probabilities A(n) and B(n) 
are independent of the waiting times. The continuum limit of this equation is obtained by using 
the replacement p^n) -> Pi (x), where Pi(x) dx is the probability of finding the particle after the 
ith jump in the interval (x, x + dx). Similarly, A(n) -> A(x) and B(n) -> B(x), with the normalisation 
A(x) + B(x) — 1. In addition we have A(n + 1) -> A(x + a) and B(n + 1) -> B(x + a). We now 
expand Eq. (156) in a Taylor series in a, a typical term being 



4(n - \) Pi (n I) -> > - .f | ,K.vl/>,(.vl |( </) - ^\ | ,l(.vl/.,(.v) |'!,' - •■■ . (157) 



where higher-order terms proportional to a 3 , a 4 , etc., are omitted. 

The system is supposed to be close to thermal equilibrium defined by a temperature T, 
A(x) ~ B(x) ~ 1/2, and to be according to detailed balance A(x) — B(x) ~ aF(x)/(2k B T). In this case 
we obtain from Eqs. (157)— (158) in the continuum limit 



Pi+i(x) = Pi (x) + — 



dx 2MX) dxk B r 



+ ■■■ . (158) 
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Rewriting Eq. (155) leads to 

.... , 1 — w(w) 1 — w(u) " . 

u u i r J 1 



(159) 



where the continuum approximation W„(u) -> W(x,u) has been made. Inserting Eq. (158) into 
Eq. (159), and using p (x) = d(x), one finds 

W(x,u) = — ^d(x) 



+ 



1 — w(w) 



!= 1 



a 2 a 2 



2 8x 



Fjx) 
/c R T 



+ 



Notice that according to Eq. (155) 
^ vv(m) 00 

X Pi-iC^M")' = W(x,u)w(u) 

u i=1 

and hence from Eq. (160) one obtains 
^, M ) = 1 - W( V ) 



w(u) f . (160) 



(161) 



+ 'Muy{ W(x, ") + y ") - y ^ 



W(x,u) 



Fjx)" 
/c R T 



+ 



(162) 



Introducing the waiting time probability density function, which for small m behaves as 

w(u) = 1 - {mf + Ci (ut) 2x + ••• (163) 

with < a < 1, the first moment of the waiting times diverges. Inserting Eq. (163) into Eq. (162) one 
is led to 

(mf - c,(wr) 2a + ••• 
W(x,u) = — — S(x) + [1 -(mf + Cl (m) 2x + •••] 



x < W(x, ") + y W(x, u) - y 



W(x, u) 



Fjx) 
k B T 



+ 



(164) 



Consider the limit a -> 0. In the standard diffusion approximation (i.e., a = 1) such a limit is 
meaningful only when both the mean waiting time and the lattice spacing a approach zero. For 
those cases where the mean waiting time diverges, the standard limit of the diffusion approximation 
breaks down. We take a -> and t -> 0, while the ratio 



2 Km ^ = 



(165) 
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is kept finite. When a = 1, K t = a 2 /(2<r» and <t> = t is the finite mean waiting time, as expected 
for this normal case. Multiplying Eq. (164) by t^w 1 ~* and using the limiting procedure defined in 
Eq. (165), we find that 

u W(x, u) - 8(x) = u 1 ~ *L FP W(x, u) , ( 1 66) 

where 

L FP =K^--— j (167) 

is the well known Fokker-Planck operator. Eq. (166) can be rewritten in t space in terms of the 
fractional Riemann-Liouville operator [232] as 

dW 

— = Dr*L FP W(x,t) . (168) 



5. 7. A fractional Fokker-Planck equation for Levy flights 



The FFPE (150) for pi = 2, i.e., for the case with finite jump length variance I 2 , reduces to the 
subdiffusive FFPE (101) discussed in Section 5.2. Here we briefly discuss the opposite case of finite 
characteristic waiting time T, but diverging I 2 -> co . This Markovian analogue of the FFPE (150), 



QW 



"8 V'(x) 
Qx mf]i 



W(x, t) 



(169) 



describes Levy flights in the force field F(x). It can be independently derived from the Langevin 
equation [159,160,216,217] 



d F(x) . . 

—x (t) r(t) 
dt mrji 



(170) 



with the Levy noise r(t) the distribution of which, p(F), is given through p(k) = exp( — K^kl") in 
Fourier space. The FFPE was investigated in some detail by Jespersen et al. [217]. 

The basic feature of the FFPE (169) is that it describes systems far off thermal Boltzmann 
equilibrium. Thus, for a harmonically bound particle underlying the potential V(x) = imco 2 x 2 , the 
stationary solution 



WJx) = 



1 



H2I2 



{x^fico^ 



(1, 1), (l,/i/2) 
(1 ?(U ), (1,^/2) J 



is Levy stable, its Fourier transform being 

n 1 mK>i\k\'' 



(171) 



WJk) = exp 



/ico 



(172) 
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Consequently, the asymptotic behaviour 
K\n Y m 



WJx) 



fico 2 \x\ 



(173) 



is enforced for large |x|. Especially, the mean squared displacement diverges, even for the stationary 
case. 

Making use of the method of characteristics, it can be shown that the solution of the FFPE (169) 
for the force F(x) = — co 2 x + F is given by 



W m j = W n [ x - 



CO 



1 — exp 



co 2 t 



flCO 



1 — exp 



co 2 t 



where W (x, t) denotes the force-free solution 



W (x,t) = — H\$ 



1 

H\x\ 



_(K l it) llfl 

which is again Levy stable [217]. 



(UM(U/2) 

(1,1), (1,1/2) J 



(174) 



(175) 



5.8. A generalised Kramers-Moyal expansion 

Taking into account higher-order terms in the Taylor expansions of the type (139), leading to the 
FPE (94), the Kramers-Moyal (KM) expansion 



QW 



n= 1 



8 

8x 



D {n \x)W(x,t) 



can be obtained, where the KM-coefficients are defined through 

f (n) w^^Mx) + (-irB(x)] . 

Alternatively, the KM-expansion foots on an expansion of the distribution function [36] 



P(x, t + x\x, t) 



dyd(x-y)P(y,t + T\x',t) , 



(176) 



(177) 



(178) 



P denoting the transition probability from x' to x during the time span x, in combination with the 
formal expansion 



d(x-y)= X 



n = 



(y-x')f _Q_ 
8x 



8{x' - x) 



(179) 



Note that in the full KM-expansion (176) no limits have to be taken as the full Taylor expansion 
is included. This is connected with the Pawula theorem [36] as either the Taylor expansion has to 
be terminated after the second order, or no proper limit can be defined. A truncation of the 
KM-expansion after the nth term, n > 2, may lead to negative solutions for the pdf W(x, t) [36]. In 
the random walk derivation, the lack of appropriate limits is obvious, as the limit is only properly 
defined for the quotient (Ax) 2 /Af. All higher terms are, however, of order (Ax) 2 + n /At. 
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A generalised KM-expansion has been obtained recently for systems underlying Levy jump 
statistics with index 1 < \x < 2 [214]: 



m - D 1 * 



00 oo p. 

X . aa D"/D^\x) + X ( - 1)» DT^"\x) 

_n=l n=l 0X 



?) \ 1+2n 

+ E o (-i) 1+n ^) &»\x) 



W(x,t) (180) 



where the generalised KM-coefficients are defined via 

DW( x )=iW- 1 — , (181) 

. / nun 

D<»>(x) = i [ "" ] (l - 5 2 J— —lA(x) - B(x)] ^ '- ^ , (182) 



2n!sin( -[1 + n\i\ 



1 + In 

ocl 



D<»»(x) = ^[#) - B(x)] ^ ^ - (183) 

Note that the occurrence of possibly imaginary coefficients in higher-order terms in the genera- 
lised Kramers-Moyal expansion (180) is due to the differentiation theorem of the Fourier 
transformation (compare footnote 16); thus, the Riesz-Weyl operator occurs to orders nfi, 
n = 1,2,... . 

A crucial feature of the generalised KM-expansion is that the lowest-order contribution 
involves a first-order spatial derivative 8/8x, and thus preserves the physical drift character. That 
means that, to lowest order, the external force in the generalised KM-expansion leads to a transla- 
tion of the pdf W(x, t). The latter statement is different from the findings of Zaslavsky et al. 
[187,193,211] who assume a generalisation of relation (179) in their description of chaotic 
Hamiltonian systems. 



6. From the Langevin equation to fractional diffusion: microscopic foundation of 
dispersive transport close to thermal equilibrium 

In this final section we briefly review a physical scenario giving some insight into the origin of the 
fractional Fokker-Planck equation for multiple trapping systems. From the continuous time 
version of the Chapman-Kolmogorov equation combined with the Markovian Langevin equation 
of a damped particle in an external force field, a fractional Klein-Kramers equation is derived 
whose velocity averaged high-friction limit reproduces the fractional Fokker-Planck equation, and 
explains the occurrence of the generalised transport coefficients K a and n a . 
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6.1. Langevin dynamics and the three stages to subdiffusion 

In his treatment of the Brownian motion of a scalar test particle in a bath of smaller atoms or 
molecules exerting random collisions upon that particle, Langevin [64] amended Newton's law of 
motion with a fluctuating force. On the basis of the resulting, stochastic, Langevin equation, the 
corresponding phase space dynamics is governed by the deterministic Klein-Kramers equation 
[35,36,81,82,278,279]. Its solution, the pdf W(x,v,t) to find the test particle at the position 
x, ... ,x + dx with the velocity v, ... ,v + dv, at time t, describes the macroscopic dynamics of the 
system. Thereby, two limiting cases can be distinguished, these being the Rayleigh equation 
controlling the velocity distribution W(v, t) in the force-free limit, and the Fokker-Planck equation 
from which the pdf W(x, t) can be derived. 

Fractional Fokker-Planck and Klein-Kramers equations have been derived and discussed for 
Levy flights which are Markovian but possess a diverging mean squared displacement. Although 
the fractional Fokker-Planck equation can be derived from the generalised master equation or 
continuous time random walk models as shown in the preceding section, a foundation on 
microscopic dynamics within the Langevin picture sheds some light on the coming into existence of 
fractional dynamics as is briefly shown in this section. 

The three stages of this model comprise the following steps: Firstly, the Newtonian motion of the 
scalar test particle experiencing a random force, in accordance to the Langevin equation 



m 



<Px 
dt 2 



= — mnv + F(x) + mr(t), v = 



dx 
dt 



(184) 



with the ^-correlated Gaussian noise r(t). Secondly, its combination with kinetic energy-con- 
serving trapping events which are ruled by the broad waiting time statistics according to 
w(t) ~ x x /t 1+x . During a trapping event the particle is temporarily immobilised. And, thirdly, 
the macroscopic average in which the long-tailed trapping events win out in the competition 
with the Langevin motion events of average duration t*, in the spirit of the generalised 
central limit theorem. This model offers some physical insight into the origin of fractional 
dynamics for systems which exhibit multiple trapping such as the charge carrier transport 
in amorphous semiconductors [95-97], or the phase space dynamics of chaotic Hamiltonian 
systems [280]. 

After straightforward calculations basing on the continuous time version of the Chapman- 
Kolmogorov equation [164,219] which are valid in the long-time limit t P max{r, t*}, one obtains 
the fractional Klein-Kramers equation [219,220] 



QW 



n 1 ~* 



J 8 

* 1 

9x dv 



n*v 



F*(x)\ fc B T 9 2 
1 + rf 



m 



m dv 2 



W(x, v, t) 



(185) 



Hereby, the Klein-Kramers operator in the square brackets has the same structure as in the 
Brownian case, except for the occurrence of the starred quantities which are defined through 
v* = v3, n* = t]3, and F*(x) = F(x)3 whereby the factor 3 is the ratio 3 = t*/t* of the inter- 
trapping time scale t* and the internal waiting time scale x. This rescaling automatically reveals the 
generalised Einstein relation K a = k B /mrj x [219,220]. 
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6.2. The fractional Klein-Kramers equation and related transport equations 



Integration of Eq. (185) over velocity, and of v times Eq. (185) over velocity results in two 
equations whose combination leads to the fractional equation [219,281] 

dW 

8x 2 



dt 



jx mn x 



W(x, t) . 



(186) 



Eq. (186) is of the generalised Cattaneo equation type [283-286] and reduces to the telegrapher's 
type equation found in the Brownian limit a = 1 [281]. In the usual high-friction or long-time limit, 
one recovers the fractional Fokker-Planck equation (101). 

The integration of the fractional Klein-Kramers equation (185) over the position coordinate, 
leads in the force-free limit to the fractional Rayleigh equation [282] 



QW 



*n* 



~_8_ k B T a 2 
Qv m Qv 2 



W(v,t) 



(187) 



Its solution, the pdf W(v, t), describes the equilibration of the velocity distribution towards the 
Maxwell distribution 



WJv) 




(188) 



The presented model for subdiffusion in the external force field F(x) = — V'(x) provides a basis 
for fractional kinetic equations, starting from Langevin dynamics which is combined with long- 
tailed trapping events possessing a diverging characteristic waiting time T. This combined process 
is governed by the long-tailed form of the waiting time pdf, manifested in the fractional nature of 
the associated Eq. (185). Physically, this causes the rescaling of the fundamental quantity r\ by the 
scaling factor 3, to result in the generalised friction constant 



t]a = rj/3 . 



(189) 



It is interesting to note that a similar process depicting a force-free trapping-walk scenario on 
a kinematics level was described in Ref. [287] in (x, ^-coordinates, revealing the subdiffusive mean 
squared displacement <x 2 (t)> oc t*. 

The Langevin picture rules the Markov motion parts in between successive trapping states. On 
this stage the test particle consequently obeys to Newton's law, in the noise-averaged sense defined 
above. Conversely, averaging the fractional Klein-Kramers equation (185) over velocity and 
position coordinates, one recovers the memory relation (d/dt)<^x(t)^> = 3 D^ x ^v(t)y between the 
mean position <C*(0^ an d the mean velocity C u (0^- This "violation" is only due to the additional 
waiting time averaging which camouflages the Langevin-dominated motion events. 



7. Conclusions 

Roughly a hundred years have elapsed since the advent of random walk and diffusion theory. 
The success of the framework developed by its most important contributor, Albert Einstein, 
obtained an additional thrust when the experimentalist Perrin came up so successfully with his 
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determination of the Avogadro-Loschmidt number. The ideas of how random walks can be used as 
a model for the transport dynamics in physical systems are still the same, and have become a joint 
venture of mathematicians, physicists, chemists, engineers, earth and life scientists. The extension of 
random walk theory to incorporate generalised statistics which no longer follow the central limit 
theorem, and memory effects violating the Markovian nature of early days random walks, has 
created a very rich tool, rich enough to be able to describe certain features of complex systems. 

Random walks are a very convenient tool in drawing a physical picture of the processes 
underlying the dynamics of systems which are of a probabilistic nature. To deal with problems 
involving boundary values or external fields, the differential equation approach is more convenient, 
however. In this report we have demonstrated the versatility of fractional equations of the diffusion, 
diffusion-advection, and Fokker-Planck type, generalising their standard counterparts. Appropri- 
ate methods of solutions such as the Fourier-Laplace technique, the method of images, and the 
method of separation of variables have been discussed. Mellin transformation techniques, term- 
by-term inversion, as well as the method of characteristics have been discussed elsewhere. We 
presented an integral transformation between the Brownian solution and its fractional counterpart. 
Moreover, a phase space model was discussed which explains the genesis of fractional dynamics in 
trapping systems. These issues make the fractional equation approach a powerful instrument. 

Note that whereas we showed the solution of the FDE for simple boundary values, the opposite 
problem, a standard equation with fractal boundary conditions, has also received some interest 
[288-293]. With the fractional approach it should be possible to discuss such type of boundary 
value problems also for anomalous diffusion. 

The main emphasis in the present work has been laid on subdiffusion. A basic feature arising in 
that context is the replacement of the exponential decay of modes by the Mittag-Leffler pattern. 
This feature involves the slow decay of the initial condition, slower dispersion, memory effects, and 
consequently a relatively slow approaching of the stationary state. Such strange dynamics have 
recently been discovered in protein systems [294]. In fact, it is supposed that dispersive kinetics are 
responsible for relaxation processes in protein systems [295,296]. 

A further advantage of the fractional approach is the relatively straightforward way of calculat- 
ing the moments. Being very simple for the CTRW approach in force-free diffusion, this advantage 
is obvious in such cases where a non-linear external force acts upon the test particle. Thereby, 
integration over the underlying fractional equation produces an ordinary fractional differential 
equation from which the moments can be inferred. 

It has been shown that subdiffusive phenomena described by the FFPE are close to thermal 
equilibrium which is reached via the aforementioned Mittag-Leffler pattern. The calculation of the 
generalised Einstein-Stokes-Smoluchowski relation, and the validity of the second Einstein 
relation in a way justify the use of the subdiffusive FFPE as direct generalisation of the Brownian 
analogue. Both phenomena have been corroborated experimentally for subdiffusive systems. 



Note added in proof 

Recently, Barkai and Silbey have presented the fractional Klein-Kramers equation [66] 
dW(x,v,t) W F(x)W nl XT 

£T + v—+ r-=f/ a oA *L FP W(x,v,t) 190 

or Ox m ov 
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for < a < 1, with the Fokker-Planck operator 

9 K B T Q 2 
L F J> = + —2 

ov m ov 

and \Y\ a ~] = s~*. This equation corresponds to the force-free mean squared displacement 



and thus describes the transition from ballistic to sub-ballistic superdiffusion. Further studies of 
this equation, especially on the non-negativity of the pdf are under way. 

We would further like to point out that a fractional Fokker-Planck equation was proposed by 
Jumarie [231], and that non-integer order spatial derivatives were presented by Onuki [246] and 
Suzuki [306]. A Levy approach to quantum pdfs has been reported by West [307]. A first passage 
time study in anomalous diffusion on the basis of the FFPE has been accomplished by Rangarajan 
and Ding [308]; we thank these authors for sending us a set of unpublished manuscripts. 
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Appendix A. Fractional differentiation and integration 

As early as in 1695, when Leibniz and Newton had just been establishing standard cal culus, 
Leibniz wrote in a letter to de l'Hospital: "Thus it follows that d 1/2 x will be equal to x Ijdx'.x 
an apparent paradox, from which one day useful consequences will be drawn." In the course of 
time, many famous mathematicians worked on this and related questions, creating the field 
which is known as fractional calculus today. The list of prominent names includes Laplace, 
Lacroix, Cayley, Pincherle, Holmgren, Grunwald, Krug, Heaviside, Laurent, Riemann, Liouville, 
Hardy, Hadamard, Levy, Weyl, Erdelyi and Kober, thus indicating the relevance of the topic in the 
mathematicians' point of view. 




(192) 
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Throughout this report, we employ the so-called Riemann-Liouville fractional integral defined 
through [226,232-235] 



t oD t - x f(t) = —\ dt' u J \' , (A.1) 



a direct extension of Cauchy's multiple integral for arbitrary complex a with Re(a) > 0. A fractional 
derivative is then established via a fractional integration and successive ordinary differentiation 
according to 

t Mf(t) = ^ nt M- n f(t) (A.2) 

with Re(/?) > 0, and the natural number n satisfies the inequality n > Re(/?) > n — 1. Through this 
definition it is clear that the fractional differentiation is non-local. As in standard calculus, 
integration and differentiation are not commutative, i.e., 



_d_ 
dx 



dt/(t) =f(t), 



X dtj t f(t) =f(x) + c , (A.3) 



certain composition rules for successive fractional operations have to be formulated, see below. 

In respect to the multiple Cauchy integral, a fractional integral may be viewed as a non-integer 
multiple integral, i.e., an a-fold integral with a a real, or even complex number. This might seem an 
odd statement; however, the same paradox is encountered in the notion of a fractal dimension 
which is a generally accepted mathematical notion. Indeed, there exist tight relations between 
fractional differintegrals and fractal geometry dimensions, as can be anticipated by the relation 
(A. 16) established below, which demonstrates that D* is a map relating power-laws of different 
indices. Fractional differintegrals even serve as sensors for the fractal graph dimensions of certain 
functions [297,298]. 

Two special cases are of importance in this report, these being the Riemann-Liouville operator 
D( for t =0, and the Weyl operator - m D% for t = — go in Eq. (A.l). Mathematically, 
the expression D?f(t) = (l/r(a))f ~ 1 *f(t) is a Laplace convolution, whereas - m D%f(x) = 
(l/r( ( u))x' 1 ~ 1 */'(x) represents a Fourier convolution. Therefore, Laplace and Fourier transforma- 
tions will be a useful tool in solving fractional order differential and integrodifferential equations. 

A.l. The Riemann-Liouville fractional operator 

The Riemann-Liouville operator, £f, defined through 



1 d" 

oD?f(t) 



fit') 
o"(f-0 



r(n - p) dt n 

for n > Re(p) > n — 1, obeys the following theorem for Laplace transformation: 

J?{ D t -«f(t)} = u-«f(u) , (A.5) 
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q > 0, which is a direct generalisation of the integral theorem for integer-order integrals. Converse- 
ly, the differentiation theorem is modified to the form [232] 

&{ Dff(t)} = u"f(u) - £ u>cj (A.6) 

which involves the quasi initial value terms 

Cj = lim oDr'-fit)- (A.7) 

t->o 

Finally, the Riemann-Liouville operator obeys the following composition rules [232]: 

D? D?f(t) = D? + Q f(t) if Q < v Q < 1 and/(0) bounded , (A.8) 

D? D?f(t) = oD? +Q [f(t) - t Cjt Q ~ k )> Q > (A.9) 

with 

^I^^TTi)- (A - 10) 

y4.2. TAe Riesz/Weyl fractional operator 

The Weyl fractional operator -^Dx has a simpler behaviour under transformations, as due to 
t = — co, no initial values come into play. Thus, its Fourier transform is 

^{_ C0 Z)£/(x)}=(ifcy/(fe) > (A.11) 

and the composition rules become 

«M «>Dl = - m Di + v \/fi,v. (A.12) 



Note that, in the main text, we preferred the simpler notation suppressing the imaginary unit, 
defined through 

^{-nDifix)} = -\k\"f(k) (A.13) 

which has somehow established in fractional applications, for instance compare [192]. In one 
dimension, the Weyl operator is equivalent to the Riesz operator which preserves the property 
(A.13) to higher dimensions. We refer to the symbol -^D^ as the Riesz/Weyl operator. 

A. 3. Differintegrable functions and an equivalent definition 

Following Oldham and Spanier [232], we define the class of differintegrable functions as all 
functions f(t) which can be expanded as a differintegrable series according to 

00 

f{t)={t-xyY.m-^ ln (A- 14 ) 

j = o 
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where t ^ 0, n is a natural number, and p > — 1. Most of the special functions used in mathemat- 
ical physics are subsumed under this definition. The definition (A. 14) of differintegrable functions 
actually is very strict, and many functions or distributions like the Heaviside jump distribution or 
the Dirac delta distribution can be differintegrated to fractional order. 

A definition equivalent to the Riemann-Liouville fractional integral (A.l) for Re(a) > is given 
through the limit [232] 



t0 D t -«f(t)= lim 



JV-i-oo 



t - t, 

~N 



o / ( / - y.) f ^ -ji +.//« 



N 



(A.15) 



Differently to Eq. (A.l), this definition also holds for fractional differentiation, i.e., Re(a) < 0. Again, 
the non-local character becomes clear. 



A. 4. Examples 

The fractional Riemann-Liouville differintegration of an arbitrary power for t = is given by 



D y t t> 



r(i+n- v) 

which coincides with the heuristic generalisation of the standard differentiation 
d"t m ml 



dt" (m-n)\ 



(A.16) 



(A.17) 



by introduction of the Gamma function. An interesting consequence of Eq. (A.16) is the non- 
vanishing fractional differintegration of a constant: 



oA v l 



1 



m - v) 

The Riemann-Liouville differentiation of the exponential function leads to 

t" 



0J D ( V = 



1^(1,1 -v,0 



(A.18) 



(A.19) 



r(i - v) 

involving the confluent hypergeometric function 1 F 1 [236]. This result can be found easily by 
differentiating term by term in the exponential series according to Eq. (A.16). On the other hand, 
for the Weyl fractional operator _ m DJ, the fundamental property of the exponential function, i.e. 
(e ( ) (n) = e' carries over to fractional orders: 



A. 5. The singular nature of the fractional operator 



(A.20) 



Some caution has to be paid to the singularity at t = t'. Consider the derivation of Eq. (A.16). 
With the definition 



B(b + l,d+ \)x 



b + d + l 



r(b + i)r(d + l) 
r(b + d + 2) : 



.b + d + l 



dt(x - tft 



b+d 



dtt b (x-t) d (A.21) 



R. Metzler, J. Klafter / Physics Reports 339 (2000) 1-77 



61 



of Euler's Beta function B, the Riemann-Liouville differintegration of the power-law can be 
achieved. The definition of the Beta function is valid for b, d > — 1 only. The fractional differenti- 
ation of order 1 — y, < y < 1, of the power x p is then 



Dl y x p 



dt(x -ty 1+y t p 



d_J_ 

dirty) 



B(y,p + l)x 1 - y+p 



Ap + l) P -i + J 



r(p + y) 

It is tempting to evaluate the parametric differentiation d/dx explicitly, according to 

df(x,t) 



_d_ 
dx 



dtf(x,t) =/(x,x) + 



dt- . 

o 9^ 



(A.22) 



(A.23) 



This leads to a term with a pole and a second term containing an integral which violates the 
definition of the Beta function. On closer inspection, this integral itself has a pole which in some 
way compensates the pole of the first term. In fact, a similar situation arises for Abel's integral 
equation: 



x/(x) = 



dt(x-ty" 2 f(t) 



(A.24) 



solved by/(x) = x~ 3l2 e~ n/x . Whereas the differentiation of the left-hand side is no problem, on the 
right two singularities "compensate" each other. The definition of fractional differintegrals 
thus involves singular integrals. However, the mathematical framework is well defined if caution 
according to above considerations is paid. 



Appendix B. Special functions: Mittag-Leffler and Fox functions 

B.l. The Mittag-Leffler function 

The Mittag-Leffler function [299,300] is the natural generalisation of the exponential function. 
Being a special case of the Fox function introduced below, it is defined through the inverse Laplace 
transform 

E,(-(,/r)-) = ^-{ u + i ! v .. }, (B.1) 
from which the series expansion 

£.( " (tfrf) - I ( m ( l /Tr) " (B.2) 
„to r(l + an) 

can be deduced. The asymptotic behaviour is 

^-(^-((t/Tmi-c 1 (b.3) 
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log 10 E 1/2 (-[t/T]^) 




Fig. 22. Mittag-Leffler relaxation. The full line represents the Mittag-Leffler function for index 1/2. The dashed lines 
demonstrate the initial stretched exponential behaviour and the final inverse power-law pattern. 



for t > t, < a < 1. Special cases of the Mittag-Leffler function are the exponential function 

£i( - t/x) = Q- tlx (B.4) 

and the product of the exponential and the complementary error function 

£ 1/2 (-aA) 1/2 ) = e" r erfc(a/T) 1 / 2 ). (B.5) 

We note in passing that the Mittag-Leffler function is the solution of the fractional relaxation 
equation [20,21,180] 



dt 



= - T\Dr a d>{t) . (B.6) 



As already mentioned in the main text, and by Glockle and Nonnenmacher [21], the Mittag- 
Leffler function interpolates between the initial stretched exponential form 

E. ( - ((/t) .,.exp(- 7 |^) (B.7) 

and the long-time inverse power-law behaviour (B.3). The Mittag-Leffler function E ll2 ( — (t/x) 112 ) 
is displayed in Fig. 22. 

B.2. The Fox function 

The Fox function, also referred to as Fox's //-function, //-function, generalised Mellin-Barnes 
function, or generalised Meijer's G-function, has been applied in statistics before it was introduced 
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into physics by Schneider [155,174,175,301] as analytic representations for Levy distributions in 
x-space, and as solutions of fractional equations. The importance of the Fox function lies in the fact 
that it includes nearly all special functions occurring in applied mathematics and statistics, as its 
special cases. Even sophisticated functions like Wright's generalised Bessel functions, Meijer's 
G-function or Maitland's generalised hypergeometric function are embraced by the class of 
Fox function. 

In 1961 Fox defined the //-function in his studies of symmetrical Fourier kernels as the 
Mellin-Barnes type path integral [237,240-243] 19 



H™(z) = H 



(a P ,A P )~ 
(b q ,B q )_ 
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dsx(s)z s 



(B.8) 



with the integral density: 

Uinbj - BrfUl W-aj + Ajs) 



ni + iHl - bj + Bjs)Y]> +1 r(aj - A jS ) 



(B.9) 



Note that the path integral in Eq. (B.8) represents just the inverse Mellin transform of %(s) [237]. 

Due to the structure of the defining integral kernel x(s) from Eq. (B.9), the Fox functions fulfils 
several convenient properties, three of which we list here: 

Proposition B.l. For k > 
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Proposition B.2. 
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Proposition B.3. The fractional differintegral of the Fox function is a map into the function class 
[180,237]: 



a Tjm,n 



(azf 



(a p ,A p )~ 
(b q ,B q )_ 



= z 



a — v jjm,n + 1 



iip+i 



q+l 



(azY 



( - a,0),(a p ,A p ) ' 
(b q ,B q ),(v-a,P). 



(B.12) 



19 According to Ref. [243], these functions have been known since at last 1868, and Fox rediscovered them in his 
studies. Their introduction into physics is due to Bernasconi et al. [155] in the study of conductivity in disordered 
systems. Schneider [301] demonstrated that Levy stable densities can be expressed analytically in terms of Fox functions. 
Wyss [174], and Schneider and Wyss [175] use Fox functions for the solution of the fractional diffusion equation. 
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An H-function can be expressed as a computable series in the form [237,242] 



U7=uj*H r (bj ~ Bj(b h + v)/B h ) m =1 r(l -a j+ Aj(b h + v)/B h ) 



H mn (z) = y y 

pq{ ' „t\ v % ni^ + i^l - bj + Bj(b h + v)/B h ) UU + i r (dj - Aj(b h + v)/B h ) 
( _ iy z (b„+v)/B„ 



v\B„ 



H™(z) ~ X res( Z ( S )z s ) 



v = 



to be taken at the points s = (cij — 1 — v)/Aj, for j = 1, ... ,n. 
Some special functions and their Fox function representation: 



(b,l) 



1 



(i + zy 



(l-r, 1)- 
(0,1) . 



az' 



05/a,l). 



1 + az* 

Maitland's generalised hypergeometric or Wright's function: 
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_(b l ,B 1 ),...,(b q ,B q ) 



— z 



= H P ',l+ i 
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Generalised Mittag-Leffler function (£ a ,i(z) = E a (z)): 



-z) = H\i 



(0, 1) 



- z 



(0,1), (1-/5, a) J jf o r(p + <xj) 



(B.13) 



which is an alternating series and thus shows slow convergence. For large argument |z| -> oo , the 
Fox functions can be expanded as a series over the residues [241] 



(B.14) 



(B.15) 
(B.16) 
(B.17) 



(B.18) 



(B.19) 



Appendix C. Some remarks on Levy distributions and their exact representation in terms of 
Fox functions 



The fundamental importance of the normal distribution is due to the Central Limit Theorem 
which, within the history of probability theory, is a consequence of the inequality of Bienayme, the 
theorems of Bernoulli and de Moivre- Laplace, and the law of large numbers. These concepts were 
extended by the works of Paul Levy, after which the generalised normal distributions are named, 
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A.Ya. Khintchine, A.M. Kolmogorov and B.V. Gnedenko, among others. The classical books 
dealing with Levy distributions and the Levy-Gnedenko generalised central limit theorem are 
listed in Refs. [56-58,88]. 

According to Levy [57], a distribution F is stable iff for the two positive constants c l and 
c 2 there exists a positive constant c such that X given by 

c 1 X 1 + c 2 X 2 = cX (Cl) 

is a random variable following the same distribution F as the independent, identically distributed 
(iid) random variables X t and X 2 . Alternatively, if 



cp(z) = <e^> 



e lXz dF(X) (C.2) 



denotes the characteristic function of the distribution F, then F is stable iff 

(p(c 1 z)(p(c 2 z) = q>(cz) . (C.3) 

A more general definition is given by Feller [88]. Let X,X 1 ,X 2 , ...,-X„ be iid random variables 
with a common distribution F. Then F is called stable iff there exist constants c„ > and y„ 
such that 

F„=X*i= c n X + y n (C.4) 

i 

where = indicates that the random variables of both sides follow the same distribution, F. 20 
Consequently, the characteristic function according to definition (C.3) fulfils the functional relation 

(p"(z) = <p{c n z)^ , (C.5) 

which can be solved exactly. The result is 

Proposition Cl. 

<A(z) = log (p(z) = iyz - c|z|* jl + ip^ co(z, a) j , (C.6) 

where a, /?, y, c are constants (y is any real number, 0<a<2, — 1 < ft < 1, and c > 0), and 
( na 

tan— n a ^ 1 , 

co(z, a) = 1 2 (C7) 

- log|z| if a = 1 . 
[n 

a is called the Levy index or characteristic exponent. From Eq. (C.6) it can be shown that the 
normalisation factor c„ in Eq. (C.3) is w 1/a . The limiting case a = 2 corresponds to the Gaussian 



As remarked by Takayasu [91], this property may be viewed as some kind of self-similarity. 
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L 




a 



Fig. 23. Parameter space (a,/?) for stable laws. All pairs of indices inside and on the edge of the diamond shape refer to 
proper stable laws. The double line denotes one-sided stable laws (OS). The letters represent the normal or Gaussian law 
(AT), the Holtsmark distribution (H), the Cauchy or Lorentz distribution (C), and the approximate log-normal distribu- 
tion (L) close to a w 0, see text. 



normal distribution governed by the central limit theorem. For /? = 0, the distribution is sym- 
metric, y translates the distribution, and c is a scaling factor for X. Thus, y and c are not essential 
parameters, and disregarding them, the characteristic function satisfies the following proposition. 



Proposition C.2. \cp(z)\ = e 



a ^ 1 . 



(C.8) 



Thus, one can write 




(C.9) 



with the new centring constant /? which is restricted in the following region: 




(CIO) 



The resulting allowed parameter space is portrayed in Fig. 23. 



Proposition C.3. The pdff x ^(x) is the Fourier transform of cp(z), defined by Eq. (C.9): 




(C.ll) 
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Thus, 

fx,fi( x ) = fx,-ti( — x ) 
and consequently 

/a,oW =fa,o( — X) 

is symmetric in x. 



(C.12) 
(C.13) 



Proposition C.4 (Levy-Gnedenko generalised central limit theorem). For iid random variables 
X 1 , X 2 , ■ ■ ■ let Y n — Yj = i Xi . If the distribution of Y n , with an appropriate normalisation, converges to 
some distribution F in the limit n -*■ co , F is stable. In particular, if its variance is finite, F is Gaussian, 
and obeys to the central limit theorem. 

Proposition C.5. The asymptotic behaviour of a Levy stable distribution follows the inverse power-law 



/^W~OTT^ a<2 - 



Proposition C.6. For all Levy stable laws with < \i < 2, the variance diverges: 

<X 2 > -> GO 

Proposition C.7. Conversely, the fractional moments of the absolute value of x, 

<\x\ 3 ) < oo 
exist for any < S < \i < 2. 



(C.14) 



(C.15) 



(C.16) 



(C.17) 



Proposition C.8. The analytic form of a stable law is given through the Fox function [301] 

fx — £#2 2 X 

for a > 1, and with the abbreviations e = 1/a and y = (a — P)/2a.;for a < 1, one obtains the result 

( - l,l),( - 



(-£,£),(- y, y) 



(C.18) 



Some reductions of these Fox representations for special cases of /? are discussed by Schneider 
[301]. From the known theorems of the Fox function, the series representations and the asymptotic 
behaviour can be determined, thus, for a > 1, one obtains 



f*A x ) = :Z 



1 " T(l + ns) 



sin(7tny)( — x)' 



n-l 



(C.19) 



n= 1 
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and the asymptotic behaviour 
1 °° m + na) 

faA*) ~ - E , sin(7may)|x| " 1 ~« (C.20) 



n= 1 



for > a — 2. An exception is the case /? = a — 2. Conversely, for a < 1, the expansions for large 
and small x are given by Eqs. (C.19) and (C.20), respectively, except the case |/?| = a. 21 

Proposition C.9. For a = 2, ft = 0, a«a* stable density is identical to the Gaussian normal 
distribution. 

Proposition C.10. For a = 1 a«a* /? = 0, stable density is identical to the Cauchy or Lorentz 
distribution 

= k^^W < C21 > 

Proposition Gil. I/O < a < 1 and ft = — a, pdff at - a (x) = Vx < is one-sided. For instance, 
the one-sided stable density for a = 1/2 ant/ /? = — 1/2 is given by 

/i/2,-i/2« = -V*- 3/2 e- 1/4 *. (C- 22 ) 

2a./ 71 



Proposition C.12. For «o? too sma// ant/ «o^ too large x, a5 we// as a « 0, the pdf f x ,- x (x) can be 
approximated by a log-normal distribution [302]: 

/ a ,_ a (x)xiexp(^-y(logx) 2 j . (C.23) 

For a = 3/2 and ft = 0, one recovers the Holtsmark distribution which is of some use in 
cosmology [88]. The above-mentioned special cases are included in the phase diagram, Fig. 23. 
Further examples can be found in the article of Schneider [301], and in the book of Feller [88]. 

In Fig. 24, the Gaussian normal distribution is compared to the Levy stable law/ 10 (x), the 
Cauchy distribution. 

A more recent monograph dealing with stable distributions is the book by Samorodnitsky and 
Taqqu [303]. Some additional historical remarks are to be found in the textbook by Johnson et al. 



21 This fact, as well as the striking similarity between results (C.17) and (C.18) are based on the fundamental property 
[237,242] 
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/(*) 




-4-2 2 4 



Fig. 24. Comparison of the Gaussian normal pdf 7i~ 1/2 e~* 2 (dashed) and the Cauchy pdf jt~ '(1 + x 2 ) -1 , the Levy stable 
law for Levy index [i = 1. Both are normalised to unity. Note that the concentration around zero is much more 
pronounced for the Gaussian. The insert shows the double-logarithmic plot of the same functions, pronouncing the 
slower power-law decay of the Cauchy pdf which turns over to the straight line in this presentation. 



Q F , . — D 



• 

p 

Fig. 25. A roaming Brownian particle P moves along a line, and radiates light. On the moving film F, the particle leaves 
a mark each time it passes the slit. The time intervals between single marks on the film is given by the one-sided 
distribution fi/z.-i/i^)- 



[304]. An interesting and readable summary of stable laws is given by Takayasu [91]. He also 
mentions a physical example, see Fig. 25, equivalent to the Polya problem [38,305]. It leads to the 
one-sided pdf /i/2,-i/2( T ); Eq. (C.22) for the distribution of the time span x between individual 
signals. An application might be in single molecule spectroscopy [138,139]. 



Appendix D. Abbreviations used 



pdf 

CTRW 



probability density function 
continuous time random walk 
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FDE 


fractional diffusion pnnation 

11 tlL. UU llcll LA111 LI o 1 V ' 1 1 v \J UClLlVJll 


DAE 


diffusion— ad vprtion pnnation 

LI 111 LL l) 1 V / 1 1 CI LI V L 1 V / 11 V/ V_l L I LI 11 V711 


FDAE 


fractional diffusion— advection equation 


FPE 


Fokkpr— Planck pnnation 

J. L ' LV LV L. 1 A 1 cl 1 1 C IV L. LI LI CL 1 1 V ' 1 1 


FFPE 


fractional Fokker— Planck eauation 


K M-exnansion 

IV 1*1 1 ' CL 1 1 l ) 1 \ / L± 


Kramers— Moval expansion 


FKKE 


fractional Klein— Kramers equation 




ndf of inst arriving at nosition x at timp t 

L/LJ.1 VI J Llo I Ul 1 1 V lilt, CI L UUiJl 11U11 CI L L1111L/ L 


W(x t) 

¥ f l-A*} L- 1 


ndf of bein& at x at time t 

V / VI 1 V / 1 L/Wllt, CL L ./V CL L UllLw L- 


\l/(x t) 


iumn ndf considered in the decounled form \l/(x t) = w(t\X(x\ 


W \ L ) 


waitmcr timp ndf* 

W CI 1 Llllg Lllllt JJLJ.1 


X(x) 


inmn length ndf 

1 uiii yj iviitjiii i/vii 




sticking nrohahilitv Fn f94^ 

iSLlV^Jvlll^ UlUUdUlllLy ; 1 . J 




ppnpralispd diffusion ropflfiHpnt snhdiffiision 

L,L11L1 Llll JLLI LI 111 UOlUll W Llllvlull I , J\A I / VL 111 LI OIL' 11 


K" 


generalised diffusion coefficient Lew flights 

*^V/ llvl LLllLJVCL Cllll CI Ul V / 11 w V / WlllVlwll I* J — J w V V 111^11 CO 


A(x x') 


transfer kprnpl 


vv o(a; 


initial condition = 1im ^ , 1/I/Yv t\ 


W Ax) 


stationary solution l/f/^fY^ = li'm. W(x t\ 

ijLLiLiuiia.i y ouiuiiuii ' ' st V*^/ — m — > 00 ' V*^) ^ / 


E (z) 


Mitta q— T ,effler function 

lYULLUg 1 — . C I 1 IV 1 1 11 1 1 C L 1 V 1 1 




Fox's iJ-function 




Boltzmann factor 


m 


mass of test particle 


1* 


generalised friction coefficient 
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